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I.  INTRODUCTION 

Several  types  of  space  flight  propulsion  systems  have  been  developed  over  the 
years.  These  include  chemical,  nuclear,  electric  and  solar  propulsion.  The  majority 
of  space  thrusters  to  date  have  been  chemical  rockets.  Although  the  Chinese  used 
rockets  over  800  years  ago,  true  development  of  rocket  propulsion  took  place  during 
this  century  [Ref.  1].  Chemical  thrusters  give  high  thrust-to-weight  ratios,  larger  than 
unity,  and  have  been  fully  developed  in  the  form  of  space  launch  vehicles  and 
attitude  control  thrusters.  In  contrast,  other  propulsion  systems  have  been  developed 
only  to  the  proof-of-concept  stage,  and  essentially  remain  at  this  stage  of 
development.  Nuclear  propulsion  was  studied  with  the  NERVA  (Nuclear  Engine  for 
Rocket  Vehicle  Application)  thruster  in  the  1960's,  and  abandoned  [Ref.  2:pp.  517- 
519].  Electric  propulsion  flights  during  the  1960's  included  the  U.S.  SERT-1  (Space 
Electric  Rocket  Test)  in  1964  and  the  U.S.S.R.  Yantari-1  rocket  in  1966.  Solar- 
electric  propulsion  was  demonstrated  via  the  SERT-2  rocket  in  1970,  powering  the 
electric  thruster  from  power  generated  by  solar  cells.  Further  electric  propulsion 
research  flights  in  the  1980's  included  the  U.S.  Navy's  NOVA-1  satellite  in  1981.  and 
Japan's  MS-T4  satellite,launched  from  the  Space  Shuttle.  Beyond  this,  nonchemical 
thrusters  have  only  been  used  in  auxiliary  roles,  such  as  station-keeping  and  attitude 


control  on  geosynchronous  satellites.  NASA's  Project  PATHFINDER  in  the  mid- 
1980's  proposed  the  use  of  a  megawatt-level  electric  plasma  thruster  for  a  manned 
Mars  mission.  However,  development  of  this  project  was  never  funded. 

In  comparing  the  different  propulsion  schemes,  a  primary  performance  indicator 
is  specific  impulse,  defined  as  the  ratio  of  thrust  to  the  rate  of  propellant  usage,  or 
alternately,  propellant  effective  exhaust  velocity  (ue),  divided  by  the  sea-level 
gravitational  constant,  (g0). 


mu        u 
I     =  I  =  -l      sec  (1) 

*&       g0 


Chemical  rockets  are  inherently  limited  in  performance  by  the  total  energy  available 
in  the  fuel/oxidizer  combustion  process,  so  that  the  total  enthalpy  available  for 
conversion  into  exhaust  kinetic  energy  is  limited.  Exhaust  velocity  is  also  limited  by 
material  heating  limitations  of  the  combustion  chamber  and  nozzle  throat,  and 
"frozen  flow  Losses"  (unrecoverable  energy  deposition  in  internal  modes  of  the  gas) 
[Ref.  3:pp.  4-5].  Peak  specific  impulse  for  liquid  chemical  propellants  is  presently  on 
the  order  of  450  seconds.  This  capability  is  completely  sufficient  for  the  tasks  of 
launch  to  low  earth  orbit  (LEO),  attitude  control,  station  keeping,  and  such  missions. 
However,  for  missions  such  as  manned  interplanetary  exploration,  chemical 
propulsion  can  be  shown  to  be  clearly  inadequate.  A  comparative  analysis  of  a  Mars 


mission  using  chemical  and  electric  propulsion  systems  shows  the  large  difference  in 
mass  payload  ratio  (final  mass/initial  launch  mass)  for  the  two  systems.  A  chemical 
system  using  a  high  impulse  Hohman  ellipse  trajectory  delivers  a  maximum  of 
approximately  10%  to  18%  of  the  launch  mass  to  the  Martian  surface  [Ref.  4:p.  115]. 
In  comparison,  an  electric  system  using  a  low  impulse  spiral  trajectory  could  deliver 
from  20%  to  60%  of  the  launch  mass,  depending  on  the  desired  transit  time.  Each 
mission  assumes  transit  from  low  Earth  orbit  to  Mars  orbit.  An  electric  propulsion 
system  would  still  need  a  high  thrust  propulsion  system  to  reach  the  Martian  surface 
Ref.  5:pp.  344-346].  The  large  difference  in  payload  ratio  is  due  to  the  much  larger 
exhaust  velocity  and  more  efficient  use  of  fuel  by  electric  propulsion.  Thus,  some 
form  of  electric  or  hybrid  electric  thruster  would  seem  to  be  in  order  for  such 
interplanetary  missions.  However,  due  to  the  low  thrust-to-weight  ratio  of  electric 
thrusters,  they  must  be  launched  into  orbit  by  other  means.  Their  usefulness  is 
restricted  to  space  thrusters,  not  to  launch  systems. 

With  specific  impulses  of  as  high  as  10,000  seconds,  electric  propulsion  offers 
the  performance  envelope  needed  for  manned  interplanetary  missions.  Electric 
propulsion  is  divided  into  three  types  of  thrusters:  electrostatic,  electrothermal,  and 
electromagnetic.  The  type  relevant  to  this  work  is  the  magnetoplasmadynamic  (MPD) 
thruster,  an  electromagnetic  propulsion  system  that  utilizes  the  Lorentz  force  created 
by  an  electric  current  together  with  its  induced  magnetic  field  to  propel  a  gas  that 
has  been  heated  to  the  plasma  state.  According  to  electromagnetic  theory,  a 
conductor  carrying  a  current  produces  an  induced  magnetic  force  perpendicular  to 


the  current.  The  applied  electric  field  and  its  induced  magnetic  field  interact  to 
produce  the  Lorentz  force  (F  =  JxB)  perpendicular  to  both  fields  on  the 
conductors.  This  briefly  summarizes  the  concept  behind  the  "self-field"  MPD 
accelerator  [Ref.  2:pp.  485-486].  MPD  performance  is  enhanced  by  adding  magnetic 
coils  to  the  thruster,  thus  strengthening  the  magnetic  field  and,  as  a  consequence,  the 
Lorentz  force  and  thrust.  This  thruster  is  appropriately  called  an  "applied-field"  MPD 
thruster.  MPD  thrusters  have  shown  specific  impulses  of  up  to  7,000  seconds  and 
efficiencies  as  high  as  70%  [Ref.  6:pp.  2-3].  Performance  of  MPD  thrusters  is  limited 
by  several  factors,  including  electrode  erosion,  current  spotting,  frozen  flow  losses, 
and  electrode  power  deposition.  Specifically,  anode  power  deposition  is  the  single 
largest  power  loss  mechanism  in  MPD  thrusters  operating  at  submegawatt  power 
levels  [Ref.  7].  In  the  following  work,  we  review  and  analytically  model  the  MPD 
anode,  including  the  sheath  and  anode  potential  drop. 


II.  LITERATURE  REVIEW 

Anode  losses  significantly  limit  magnetoplasmadynamic  (MPD)  thruster 
performance.  Much  effort  has  been  placed  on  characterizing  these  losses  and  on  the 
nature  of  power  deposition  in  the  anode  [Refs.  8-14].  As  much  as  80%  of  thruster 
total  power  may  end  up  being  deposited  in  the  anode  at  sub-megawatt  power  levels 
[Refs.  8,15].  This  power  deposition  together  with  current  constriction  at  the  anode 
surface  present  serious  problems  to  thruster  cooling  and  performance,  as  well  as  to 
anode  lifetime.  Before  any  practical  design  can  be  achieved,  a  more  thorough 
understanding  of  the  phenomena  at  the  anode,  particularly  the  anode  sheath,  must 
be  gained.  Studies  have  shown  that  the  anode  power  fraction  depends  on  thruster 
power,  current,  mass  flow  rate,  and  the  parameter  J2/rh  [Refs.  8,12,13,16].  It  has  also 
been  shown  that  the  anode  fall  voltage  is  inversely  proportional  to  anode  current 
density  [Refs.  13,16].  It  is  believed  that  a  better  understanding  of  the  role  of  an 
elevated  electron  temperature,  of  current  flow  dimensionality,  and  of  current 
unsteadiness  are  prerequisites  for  the  evolution  of  any  practical  MPD  thruster  design. 

Computer  codes  that  accurately  describe  observed  data  from  steady-state  MPD 
thrusters  have  been  developed  [Refs.  17-19].  However,  these  codes  do  not  adequately 
describe  observed  data  from  quasi-steady  thruster  experiments.  It  has  been  suggested 
that  the  lack  of  proper  electrode  modelling  (i.e.,  sheaths  and  fall  potentials)  in  these 


codes  may  explain  this  discrepancy  [Ref.  6].  Limited  analytical  work  has  been  done 
in  modelling  the  sheath  and  ambipolar  regions  at  the  anode,  influenced  perhaps  by 
the  difficult  set  of  coupled,  nonlinear  partial  differential  equations  involved.  Hugel 
[Ref.  12]  and  Subramaniam  [Ref.  20]  address  the  influence  of  the  sheath  region,  but 
do  not  model  the  electric  field,  temperature,  or  sheath  fall  voltage. 

Given  the  minuscule  extent  of  the  sheath  versus  thruster  anode  curvature,  the 
problem  at  first  appears  one-dimensional  in  nature.  A  one-dimensional,  collisional, 
equilibrium  solution  can  satisfactorily  reproduce  the  observed  electric  field  and 
charge  density  distributions  for  the  entire  sheath  and  ambipolar  regions  for  a  sheath 
where  the  electron  temperature  equals  that  of  the  heavy  species  [Ref.  21].  However, 
this  model  cannot  describe  any  decrease  in  current  density  away  from  the  surface,  or 
current  constriction,  at  the  anode  surface  which  might  be  necessary  in 
nonequilibrium.  A  two-dimensional  model,  developed  by  Biblarz  and  Dolson  [Ref. 
14],  represents  these  phenomena  and  predicts  the  voltage  drop  in  the  region.  It  is 
shown  that  the  sheath  must  account  for  a  majority  of  the  anode  voltage  drop,  and 
that  the  sheath  extent  must  be  greater  than  the  Debye  length  [Refs.  14,21].  Thus, 
a  combination  of  one-  and  two-dimensional  approaches  appears  to  better  describe 
sheath  behavior.  Incorporation  of  modelling  of  this  sort  may  improve  the  ability  of 
the  computer  codes  cited  above  to  properly  describe  quasi-steady  thrusters. 

Next,  a  description  of  the  anode  region  is  presented  in  order  to  delineate  some 
of  the  possible  effects  of  temperature. 


III.  ANODE  DESCRIPTION 

A.       THRUSTER  GEOMETRY  DESCRIPTION 

The  majority  of  plasma  thrusters  to  date  have  consisted  of  a  central  cathode 
rod  surrounded  by  an  annular  shell  anode,  as  shown  in  Figure  1  [Ref.  23].  The 
thruster  illustrated  is  sufficient  to  produce  needed  thrust  at  current  levels  above  one 
kiloamp.  Below  this  level,  an  external  magnetic  field  produced  by  an  annular  magnet 
is  needed  to  ensure  sufficient  Lorentz  force  on  the  plasma  propellant  to  meet  thrust 
requirements.  [Ref.  8].  As  illustrated  in  Figure  1,  the  JxB  body  force  simplifies  into 
an  axial  (jrB6)  body  force,  which  provides  direct  electromagnetic  thrust  ("blowing"), 
and  a  radial  ( -jzBQ)  body  force,  which  provides  electromagnetic  compression  of  the 
plasma  and  a  subsequent  pressure  force  along  the  cathode  surface  ("pumping").  [Ref. 

6] 

A  notable  exception  to  this  geometry  is  the  Stationary  Plasma  Thruster  (SPT), 
a  design  from  the  former  Soviet  Union.  The  SPT  is  an  example  of  a  plasma 
propulsion  system  known  as  a  Hall  Current  Plasma  accelerator.  An  electric  field  is 
applied  axially  to  a  stream  of  flowing  plasma,  in  addition  to  a  magnetic  field  with  a 
strong  radial  component,  which  is  applied  by  an  external  electromagnet.  When  the 
axial  electric  field  is  applied  and  a  current  flows  through  the  plasma,  an  azimuthal 
component  of  current  is  induced,  i.e.,  the  "Hall"  current. 
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Figure  1  -  Magnetoplasmadynamic  (MPD)  Thruster,  with  Axial  and  Radial  Forces 
on  Plasma  Indicated.  [Ref.  23] 


Thrust  is  produced  by  electrostatically  accelerating  the  ions  in  the  plasma,  as  well  as 
through  the  induced  Lorentz  force  mentioned  above.  A  strong  radial  magnetic  field 
is  applied  to  the  plasma,  whose  properties  are  controlled  to  make  the  electron 
Larmor  radius1  small  compared  to  the  mean  free  path2,  while  the  ion  Larmor  radius 
is  comparatively  large.  As  a  consequence,  electron  mobility  in  the  axial  direction  is 
greatly  reduced.  Thus,  the  electric  field  energy  is  given  mainly  to  the  ions,  producing 
axial  ion  acceleration.  Collisions  with  neutral  particles  serve  to  accelerate  the  entire 
neutral  plasma.  [Ref.  24] 

A  pair  of  the  final  prototype  design  developed,  the  SPT-100,  have  been 
acquired  by  NASA  recently  from  Fakel  Enterprise  in  Kaliningrad,  Russia,  and  are 
undergoing  performance  evaluation  at  the  Jet  Propulsion  Laboratory.  Designed  at 
the  Kurchatov  Institute  of  Atomic  Energy  (IAE)  in  Moscow,  USSR  in  the  1960's, 
smaller  versions  of  the  SPT-100  (SPT-50  and  SPT-70)  were  flown  beginning  in  19723. 
A  specific  impulse  of  1,600  seconds  and  50%  efficiency,  as  well  as  space  flights  of 
fifty  similar  thrusters  is  claimed.  The  specific  operational  characteristics  of  the 
thruster  are  not  well  understood  presently.  Bohm  diffusion  of  electrons  and  a 
phenomenon  called  "near-wall  conductivity"  have  been  proposed  to  explain  the 
thruster's  operation.  This  thruster  is  shown  in  Figure  2.  [Ref.  25] 


1  Larmor  radius  is  the  radius  of  the  helix  traversed  by  a  charged  particle  moving 
in  a  magnetic  field. 

2  Mean  free  path  is  the  distance  traveled  by  a  particle  before  making  a  collision. 

The  suffix  (i.e.,"-70")  indicates  the  characteristic  diameter  of  the  thruster  in 
millimeters. 


(95%) 


MAIN   DISCHARGE 
CHAMBER 


FLECTROMAONET 

MAGNETIC 
POLE   PIECE 


DIELECTRIC  WALLS 


POWER  SUPPLY- 


Figure  2  -  Stationary    Plasma  Thruster  [Ref.  25] 
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B.       ELEMENTARY  SHEATH  FORMULAE  DESCRIPTION 

1.        Discussion 

Voltage  losses  and  anode  power  deposition  account  for  most  of  the 
inefficiency  of  plasma  thrusters.  In  order  to  understand  these  losses,  the  anode  region 
must  be  understood  and  related  phenomena  explained  and  modelled.  As  shown  in 
Figure  3,  a  substantial  drop  in  voltage  occurs  in  a  short  distance  from  the  anode 
surface. 
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Figure  3  -  Electric  Field  Between  Two  Electrodes,  Including  "Positive  Column".  [Ref. 
27] 

The  anode  fall  region  may  be  divided  into  two  parts,  the  sheath  and  ambipolar 

regions.  The  plasma  attempts  to  adjust  itself  near  electrodes  so  as  to  shield  the  main 

body  of  the  plasma  from  the  electric  field  [Ref.  26].  The  sheath  is  the  region  closest 


to  the  anode  surface  within  which  the  ion  and  electron  number  densities  are  unequal, 
with  the  electrons  dominating  the  region.  A  high  electron  space  charge  exists  at  the 
anode  surface.  This  is  caused  by  the  anode  collecting  incoming  electrons  in 
completing  the  arc  current  with  the  cathode.  Positive  ions  are  produced  within  the 
sheath  by  electron  impact  of  neutral  gas  molecules,  and  the  ions  are  repelled  toward 
the  cathode.  At  the  cathode  end  of  the  anode  drop  region,  the  density  of  positive 
ions  is  high  enough  to  almost  neutralize  the  electron  space  charge,  thus  forming  the 
positive  column  or  core  plasma.  The  essential  positive  ion  current  is  created  in  this 
way  near  the  anode.  A  more  complete  description  of  this  process  may  be  found  in 
Cobine  [Ref.  27]  and  von  Engel  [Ref.  28].  A  fundamental  characteristic  of  plasma 
behavior  is  its  tendency  toward  electrical  neutrality.  Whenever  local  charge 
concentrations  arise  or  external  potentials  are  introduced  into  a  system,  these  are 
shielded  out  in  a  distance  known  as  the  "Debye  length"4.  This  distance  must  be  much 
smaller  than  the  system  dimension  for  the  ionized  gas  to  be  considered  a  plasma 
[Ref.  29].  Equation  (2)  gives  the  Debye  length  [Ref.  26]. 
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The  Debye  length  effectively  describes  the  radius  of  a  shell  around  a  charged 
particle  outside  of  which  the  potential  of  the  particle  is  not  seen. 
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Another  distance  of  interest  is  the  electron  mean  free  path,  or  distance  traveled  by 
a  particle  before  making  a  collision.  Equation  (3)  is  from  a  derivation  of  Lin,  Resler, 
and  Kantrowitz  [Refs.  30,31]  giving  the  mean  free  path,  with  Xs  being  the 
approximate  sheath  length. 


k  =  0.12 


1 


n(e2/3kT)2ln(Xe2/3kT) 


(3) 


Since  the  sheath  extends  at  most  a  few  mean-free  lengths  from  the  anode  surface, 
curvature  of  the  anode  does  not  affect  the  governing  equations  in  high  pressure 
discharges.  Thus,  the  region  may  be  described  in  one  dimension,  the  distance  "y" 
from  the  anode  surface.  While  the  Debye  length  is  sometimes  assumed  as  the  sheath 
extent,  Reference  22  showed  that  the  sheath  thickness  is  a  function  of  the  anode  fall 
voltage  and  the  electron  temperature.  Equation  (4)  gives  the  appropriate  form. 
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An  example  case  with  a  fall  voltage  of  100  volts  gives  a  sheath  extent  of 
Xs  =  2.352xl0"5  m.  This  compares  to  a  computed  Debye  length  of 
XD  =  l.690xl0"6  m.  Therefore,  the  sheath  can  be  an  order  of  magnitude  larger 
than  the  Debye  length. 
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Nasser  [Ref.  32]  discusses  an  elementary  theoretical  approach  to  the  glow 
discharge  problem.  He  suggests  a  set  of  four  one-dimensional  ordinary  differential 
equations,  including  the  electron  and  ion  current  and  number  density  equations,  in 
addition  to  Poisson's  equation.  Most  solution  attempts  have  failed,  with  the  boundary 
conditions  being  identified  as  the  culprit.  A  similar  attempt  for  the  plasma  thruster 
is  discussed  below. 

2.       Simplified  Formulation 

The  steady  probe  equations  are  first  written  [Ref.  21]  in  their  simplest 
form.  The  anode  is  assumed  to  operate  as  a  heavily  biased  probe,  which  is  true  for 
low  enough  currents  when  the  anode  is  not  a  source  of  ions.  Whenever  the 
temperature  can  be  considered  fixed,  the  energy  equations  are  implicitly  satisfied 
and,  since  ion  inertia  is  neglected,  the  resulting  set  consists  only  of  two  species 
continuity  equations  and  Poisson's  equation.  These  equations  are  written  in  terms 
of  y,  which  is  the  coordinate  outward  from  the  planar  positive  surface.  Constants  and 
variables  are  listed  in  Table  1. 
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TABLE  1  -  NOMENCLATURE 


a...characteristic  length  of  plasma 
D±  e... species  diffusion  coefficient 
e... elementary  charge  constant 
E... electric  field 


n.... species  number  density  at  core  plasma 

N... total  number  density 

T... temperature 

T0... neutral  species  temperature 


E„... electric  field  at  anode  surface        a  ...two-body  recombination  coefficient 
E„... electric  field  at  core  plasma  e0... permittivity  constant 


jie... species  current  density 
J. ..total  current 
k...Boltzmann's  constant 
K.. .current  parameter 
nie... species  number  density 
ne...time  rate-of-change  of  ne 


v  ...ionization  coefficient 

Hi#e... species  mobility  coefficient 

<X>a... anode  fall  potential 

X... mean-free  distance 

Xd...Debye  length 

Xa... Sheath  thickness 


Note:  Species  subscripts  denote  ions  (i)  and  electrons  (e). 
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j,  =  e/inE  -  (eD)^i  (5) 

dy 


j    =  tfinE  +  (eD) — !  (6) 

dy 


^  =  -(n   -  n)  (7) 


J  =  j,  -  1  (8) 

M     -^  (9) 

'•'       kT 

Here,  the  j's  are  species  contributions  to  the  total  current  density.  The  existence  of 
negative  charges  as  free  electrons  is  pivotal  in  the  formulation.  Next,  the  Einstein 
relation,  equation  (9),  is  introduced  to  write  the  mobilities  in  terms  of  the  diffusion 
coefficients.  We  assume  that  the  diffusion  coefficients  remain  constant  in  the 
problem. 

Equations  (5)  and  (6)  are  next  solved  for  dnle/dy.  The  species  current 
density  equations  are  found  from  the  net  reaction  rate  of  the  plasma.  Equations  (10) 
and  (11)  combine  to  produce  space  derivatives  for  species  current  density. 
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n    =  vn    -  ann  (10) 

*-  I*,  eh  (ID 

dy         dy 

Combining  equations  (5)-(ll)  produces  a  set  of  five  coupled,  non-linear  differential 
equations  describing  the  sheath.  These  are  nondimensionalized  to  adjust  all  variables 
to  the  first  order,  and  are  rewritten  below  as  equations  (12)-(16),  with 
nondimensionalized  variables  denoted  by  "x".  Nondimensionalization  can  be 
accomplished  as  follows:  The  species  number  densities  ne,n„  are  divided  by  their 
values  at  infinity  to  produce  output  from  the  anode  surface  to  unity  at  the  ambipolar 
boundary.  The  current  densities  je,j,  are  divided  by  the  total  current,  allowing  the 
output  to  show  the  "mirror  behavior"  of  the  two  currents.  The  electric  field  is  divided 
by  the  initial  anode  value  to  give  output  starting  from  unity  at  the  surface  and 
decreasing  to  the  final  core  field  value.  The  variable  "y"  is  divided  by  the 
characteristic  length5  of  the  plasma,  "a",  producing  y.  These  corrections  allow  all 
output  to  vary  in  the  range  from  zero  to  one,  as  a  function  of  distance  from  the 
anode. 


5The  characteristic  length  is  defined  so  as  to  cancel  the  multiplying  factor  in  the 
electric  field  equation,  (14),  (a  =  1.107  x  10').  This  allows  a  physical  interpretation  of 
the  ion/electron  number  densities,  as  well  as  the  decay  rate  of  the  electric  field. 
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Attempts  to  solve  this  equation  set  using  the  computer  code  discussed  below  shows 
the  set  to  be  extremely  sensitive  to  initial  conditions.  The  computer  code  solver  uses 
a  "marching"  scheme  from  the  anode  to  the  undisturbed  plasma.  The  initial 
conditions  are  chosen  to  produce  the  electric  field  potential  drop  observed  in  actual 
thrusters.  First  and  second  space  derivatives  of  the  electric  field  are  used  as 
diagnostic  checks  to  ensure  reasonable  output  values  and  indicate  instability  of  the 
integration  process.  Figure  4  shows  the  required  resulting  curves  for  the  electric  field 
and  its  first  and  second  derivatives. 
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Figure  4  -  Electric  Field  and  First  Two  Space  Derivatives  Used  as  Diagnostic  Checks 
for  Integrator  Output.  (Plotted  for  a  Generic  Exponential  Function). 
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Ecker  characterizes  the  plasma  at  the  anode  as  a  double  sheath,  with  the  inner 
section  called  the  "inertia  sheath",  and  an  outer  section  called  the  "energy  loss 
section".  The  inner  section  shows  a  potential  rise  of  the  order  of  one  volt,  with  the 
outer  section  showing  the  exponential  potential  drop  shown  in  Figure  4.  While  this 
double  sheath  may  in  fact  describe  the  actual  sheath  region,  the  formulation  above 
only  models  the  potential  drop  portion  of  the  sheath,  and  does  not  attempt  to 
produce  the  potential  rise  of  the  inner  sheath.  In  addition,  Ecker's  current 
constrictions  are  of  a  "macroscopic"  nature,  whereas  those  of  Reference  14  and  this 
work  are  "microscopic"  [Ref.  33]. 

Data  for  a  6,000°K  Nitrogen  plasma  were  used  to  test  the  equation  set  [Ref. 
21].  Producing  a  proper  solution  required  adjusting  the  initial  conditions  to  force  the 
curve  shapes  discussed  above.  Using  Equations  (2),  (3),  and  (4),  the  mean  free  path, 
Debye  length,  and  approximate  sheath  extent  are  calculated  as  X  =  9  .  3  52x10 "3  m, 
XD  =  1 .  690xl0~6  m,  and  Xs  =  2  .  352 xlO"5  m  (this  assumes  a  drop  voltage  of  100 
volts). 
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3.       Approximate  Formulation 

Reference  21  explores  the  above  equation  set  by  taking  advantage  of  the 
symmetry  among  the  equations,  and  introducing  two  parameters,  K  *  and  K " ,  shown 
below. 

K*  s  JL  +  Jf_  (17) 

eD        eD 

i  e 

-K-  =  _L   -  _L  (18) 

eD        eD 

It  can  be  shown  that  the  resulting  equations  can  be  manipulated  to  yield  a  single, 

ordinary  differential  equation  for  the  K's  in  terms  of  the  electric  field.  The  resulting 

equation  can  be  scrutinized  for  two  distinct  temperature  regimes.    Note  that  while 

the  total  current  density,  J,  is  constant  in  a  steady,  one-dimensional  case,  the  K's  can 

vary  and  will  in  turn  also  depend  on  the  degree  of  reactivity  of  the  plasma  (tie),  i.e., 

5  .  ^k  -  en  (19) 

dy        dy 

Because  ion  diffusion  is  much  slower  than  electron  diffusion,  it  can  be  shown  that  the 

K's  are  related  by 

K"  «  -K-  +  —  (20) 

eD 

e 

As  will  be  evident,  at  the  electrode  surface  the  K's  are  equal  to  each  other  and  at  the 
undisturbed  plasma,  K"  =  0.   The  total  current  density  may  be  evaluated  from 
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J  =  en  v  (21) 

00       coo  x  ' 

where  ve(S  is  the  electron  drift  velocity  beyond  the  ambipolar  region  which  is  strictly 
a  function  of  E./N,  (i.e.,  of  the  ratio  of  undisturbed  electric  field  to  the  total  number 
density). 

a.       Effects  of  Temperature  on  Anode  Constriction 

It  is  useful  to  investigate  the  overall  effects  of  temperature.  Since 
temperature  will  be  considered  constant,  it  comes  in  as  a  parameter  in  this 
formulation  whereas  charge  density  and  electric  field  remain  as  variables.  Intuitive 
arguments  will  be  introduced  which  suggest  that  the  electron  and  ion/neutral 
temperatures  play  a  rather  singular  role  in  determining  the  intrinsic  dimensionality 
of  the  problem,  (i.e.,  there  are  cases  when  the  geometry  of  the  current  lines  is  not 
necessarily  impressed  by  the  electrode  geometry).  Since  the  problem  is  described  by 
moderate  pressure,  largely  collisional  sheaths,  the  ion  and  neutral  temperatures  are 
anticipated  to  remain  reasonably  equal.  Depending  on  the  gas,  the  electron 
temperature,  on  the  other  hand,  can  be  elevated  from  the  gas  temperature  at  the 
anode  where  actual  magnitudes  depend  on  the  local  value  of  E/N.  In  order  to  get 
a  perspective  on  the  effects  of  temperature,  we  shall  consider  two  extremes,  namely, 
the  case  where  the  electron  and  ion  temperature  are  the  same  (the  equilibrium  case) 
and  the  case  where  the  electron  temperature  is  substantially  elevated  from  that  of 
the  ions/neutrals  (the  two-temperature  case). 
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(1)      Case  I:  Tt  =  T,  =  T0    (Equilibrium) 

The  charge  densities  can  be  eliminated  by  combining  equations 
(5)-(9),  (17)  and  (18).   The  resulting  equation  can  be  shown  to  be 
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(22) 


If  the  electric  field  decreases  monotonically  from  the  wall  to  the  undisturbed  plasma 
(i.e.,  from  E0  ->  EM),  then  as  y  ->  »,  E  -  E,.,     E'  -»  0,     E"  -•  0. 
So  that  in  equation  (22)  above  the  "outer  solution"  becomes: 


K*  = 


2J 
eD 


(23) 


Now  this  represents  an  acceptable  solution  from  a  physical  point  of  view.  Moreover, 


as  y 


00 


h     «  D(K7  «  0 


(24) 


which  is  also  acceptable  for  an  equilibrium  situation  at  the  undisturbed  plasma. 
Results  [Ref.  21]  are  shown  in  Figure  5  for  the  case  of  nitrogen  at  6000°K  using  an 
approximate  electric  field  distribution. 
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Figure  5  -  Electric  Field  and  Species  as  a  Function  of  y,  Distance  From  Anode.  An 
Approximation  Using  a  Shaped  Electric  Field  and  Isothermal  Plasma  [Ref.  21]. 
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(2)      Case  II:    Te  >  >  T,  =  T0    (Two-Temperature) 

In  this  case  the  same  procedure  as  before  yields  the  following 
equation  where  terms  divided  by  Te  have  been  dropped  when  compared  to  their 
counterparts  divided  by  T0. 

(K*V  -  —  E    =     2EJ     +  S  EE"  +  _fi-E"E'  -  l^E'"  (25) 

E  kT  D        kT  eE  e 

o        e  o 

Assuming  the  same  monotonic  decrease  as  before  for  the  electric  field  from  the  wall 
to  the  plasma  proper,  as  y  -*  »,  E  -»  E„,      E'  -♦  0,      E"  ■*  0. 
Then  the  outer  solution  becomes 


dK*  2eEJ  ...    .    w  A  fn^ 

with  n    >  0  (26) 


dy  kT  eD 

►  o  e 

Or,  K*  -♦   (constant)  y  +  constant,  and  he(D  keeps  increasing  with  y. 

This  is  not  the  proper  outer  solution  for  the  one-dimensional,  equilibrium 
plasma  that  we  seek  because  the  net  ionization  rate  continues  to  increase  well  inside 
the  plasma  proper  where  conditions  should  saturate,  yielding  a  constant  electric  field. 
Therefore,  as  formulated,  Case  II  is  not  amenable  to  a  one-dimensional  solution. 
References  14  and  21  show  how  this  case  can  be  analyzed  under  a  multidimensional 
approach.  These  references  also  discuss  a  method  for  describing  the  electron 
temperature  as  a  function  of  E/N,  then  how  to  couple  a  simplified  energy  relation 
which  satisfactorily  describes  a  two-temperature  plasma.  The  necessary  ingredient 
to  make  equation  (26)  approach  zero  beyond  the  decrease  of  E  to  Ex  is  to  allow  J 
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to  fan  out  as  indicated  in  Figure  6.  Thus,  in  equation  (26),  the  product  "EJ"  can 
bring  down  the  charge  production  rate  to  arbitrarily  low  values.  Alternatively,  it  is 
possible  to  explore  techniques  of  bringing  the  electron  temperature  down  to  be  in 
closer  equilibrium  with  the  ions  and  neutrals.    Transpiration  cooling  is  one  such 


means. 


U(x,y) 


B 

o 

undisturbed        plasma 


J(x,y) 


7f?  / 
1- —  node 


Figure  6  -  Two-Dimensional  Model  of  Current  Paths  Showing  Periodic  Structure. 
Thermal  Instabilities  and  Inhomogeneities  Would  Favor  One  Site  Over  Others  and 
a  Single  Macroscopic  Constriction  May  Then  Be  Produced.  [Refs.  14,  34] 
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b.       Similarity  to  Vacuum  Arc  Phenomena 

Instability  phenomena  observed  in  vacuum  arcs  [Ref.  35]  are  very 
similar  to  those  observed  in  self-field  thrusters  [Ref.  12].  After  the  establishment  of 
the  current,  the  anode  region  operates  in  a  vapor  that  issues  from  the  electrodes. 
In  vacuum  arcs,  Miller  characterizes  the  anode  region  as  operating  in  one  of  five 
distinct  modes,  ranging  from  a  passive,  low  current  mode  to  a  high  current,  fully 
developed  spot  mode  [Ref.  36].  Given  the  similarities  mentioned  above,  vacuum  arc 
anode  research  should  be  helpful  in  the  understanding  of  MPD  thruster  transition 
to  the  anode  spot  mode.  Existence  diagrams  after  Miller  [Ref.  36]  are  shown  in 
Figure  7,  which  divide  operating  modes  into  regions  as  a  function  of  anode  current 
versus  electrode  geometry.  Figure  7  shows  the  transition  from  glow  to  spot  mode. 
Anode  spot  formation  at  high  currents  is  clearly  a  factor  in  limiting  anode 
lifetime.  Various  phenomena  have  been  related  to  anode  spotting.  Hugel  [Ref.  12] 
relates  the  transition  to  spotting  mode  to  an  increase  in  J2/rh  above  a  critical  level. 
A  separate  factor  connected  with  the  spot  mode  is  surface  temperature  of  the  anode. 
Rich,  et.al,  [Ref.  37]  show  that  anode  spotting  is  preceded  by  a  luminous  "footpoint" 
and  followed  by  local  melting  prior  to  spot  formation.  Separately,  Schuocker  [Ref. 
38]  finds  a  connection  between  spotting  initiation  and  the  factors  of  anode 
evaporation  and  magnetic  constriction  in  vacuum  arcs  with  high  currents. 
Experimental  investigations  must  be  performed  to  see  if  the  above-mentioned 
vacuum  arc  criteria  apply  to  self-field  thrusters. 
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Figure  7  -  Anode  Discharge  Modes  as  a  Function  of  Current  and  Gap  Length.  [Ref 
36] 
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C.       COMPUTER  CODE 

Rather  than  using  linear  approximations  to  equations  (12)-(16),  the  nonlinear 
set  was  used,  with  initial  conditions  adjusted  in  an  attempt  to  produce  observed 
electric  fields  from  probe  data.  First  and  second  space  derivatives  of  the  electric  field 
were  used  as  diagnostic  checks  to  ensure  computed  output  was  reasonable.  Initial 
conditions  computed  from  the  approximate  formulae  in  Reference  21  were  used.  The 
equation  set  above  presents  a  difficult  problem  for  two  reasons,  nonlinearity  and 
multiple  time  constants.  The  species  number  density  equations,  (12)  and  (13),  both 
contain  a  nonlinear  term,  each  with  a  time  constant  of  its  own.  In  addition,  the 
electric  field  equation,  (14),  adds  a  possible  third  time  constant.  This  constitutes  a 
"stiff  set  of  equations.  Attempts  were  made  to  solve  the  set  with  the  data  discussed 
above,  using  Gear's  method  of  backward  differentiation,  in  hopes  that  the  variables 
would  change  slowly  enough  with  each  iteration  to  render  a  convergent  iterative 
process.  As  described  in  Reference  39,  if  some  reactions  are  slow  and  others  fast 
among  a  set  of  coupled  equations,  the  fast  ones  will  control  the  stability  of  the 
method.  This  is  addressed  in  the  DGEAR  program  available  from  the  International 
Mathematical  &  Statistical  Library  (IMSL).  The  latter  software  contains  an  Adams 
predictor-corrector  method,  as  well  as  Gear's  method,  which  is  well  known  for  its 
success  at  solving  stiff  equation  sets.  The  DGEAR  software  allows  for  a  choice  of 
functional  or  chord  iteration  methods,  as  well  as  a  choice  of  Jacobian  matrices.  A 
more  detailed  discussion  of  this  software  can  be  found  in  Reference  39  and  in  the 
IMSL  library.  [Ref.  39] 
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D.       COMPUTATIONAL  RESULTS 

Numerous  computer  runs  were  completed  using  the  initial  conditions  taken 
from  Reference  21.  In  addition,  data  for  the  ionization  coefficient  v,  Figure  8,  was 
taken  from  References  40  and  41. 
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Figure  8.  Ionization  Coefficient  v  as  a  Function  of  E/N  for  Nitrogen  for  Various 
Vibrational  Temperatures  (Refs.  40,41). 

Various  combinations  of  initial  conditions  and  ionization  coefficient  were  used.  As 

mentioned    above,    the    electric    field    and    its    space    derivatives   were    used    as 


30 


diagnostic/reasonability  checks  on  the  output.  Individual,  as  well  as  multiple 
computer  runs  were  attempted  to  model  the  sheath  region.  Nonlinearities  in  the 
equation  set  are  clearly  seen  in  Figure  9.  The  ion  number  density  does  not  reach  that 
of  the  electrons,  and  the  latter  population  growth  rate  continues  to  grow  without 
bound.  The  shape  of  the  electron  population  curve  is  very  sensitive  to  its  initial  value. 
As  shown  in  Figure  9,  the  latter  population  has  too  high  a  growth  rate  when 
compared  to  the  ion  population,  and  the  latter  does  not  "catch  up".  Increasing  the 
initial  value  of  ne  flattens  out  this  curve  to  a  reasonable  shape.  Above  an  initial 
value  of  approximately  0.06,  however,  the  plot  of  ne  "dips"  after  a  certain  distance 
and  then  continues  to  increase  as  expected.  This  gives  an  approximate  upper  value 
for  this  initial  value.  To  avoid  instabilities  like  this,  small  "slices"  were  taken  of  the 
output  after  a  small  number  of  integration  steps  and  multiple  runs  were  used  to  form 
a  "cut  and  paste"  plot  of  the  region.  When  a  reasonable  plot  shape  was  produced,  the 
value  of  ionization  coefficient  was  varied  in  the  "slices"  to  attempt  to  produce  the 
required  end  values  for  electric  field  and  species  population.  Both  multiple  and 
equilibrium  values  for  the  ionization  coefficient  were  used.  When  the  data  showed 
signs  of  instability  and  failure  to  follow  the  required  forms  of  Figure  4,  a  "slice"  was 
made  in  the  data  stream,  and  the  data  points  from  this  point  used  to  start  a  new 
computer  run.  This  approach  was  taken  in  the  hope  of  avoiding  singularities  in  the 
integration  from  anode  surface  to  ambipolar  region.  In  addition  to  the  diagnostic 
checks  shown  in  Figure  4,  an  additional  data  check  is  provided  by  the  transition  from 
the  sheath  to  the  ambipolar  region. 
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Figure  9.  Species  Number  Density  plots  for  Individual  Computer  Run,  Showing 
Divergent  Tracks  for  Ion  and  Electron  Populations,  and  Effect  of  Nonlinearity. 


As  shown  in  Figure  5,  the  species  number  densities  are  equivalent  in  this  region,  as 
are  their  change  rate.  Thus,  setting  Equations  (12)  and  (15)  equal  to  each  other  and 
solving  for  n  yields  a  value  of  0.5  in  the  ambipolar  region.  As  indicated  in  Figures 
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10-11,  the  output  produces  the  desired  plot  slopes  for  electric  field  and  species 
number  density.  However,  the  number  density  plots  cross  long  before  approaching 
the  required  value  of  0.5.  In  addition,  neither  electric  field  nor  species  number 
density  approaches  an  equilibrium  value  or  shows  sign  of  levelling  off.  Apparently, 
the  multiple  time  constants  and  nonlinear  portions  of  the  number  density  equations 
combine  to  create  a  seemingly  intractable  system.  Solutions  for  this  system  may  be 
possible  for  specific,  individual  initial  condition  sets,  but  the  problem  does  not  appear 
amenable  to  this  approach  in  general.  A  one-dimensional  system  such  as  this  may  be 
better  described  through  the  approach  of  boundary  layer  theory  or  nonlinear 
dynamics  and  chaos.  Given  the  effort  and  difficulty  involved  in  the  latter,  a  one- 
dimensional  approach  such  as  that  modelled  above  does  not  appear  useful.  A 
combination  of  one-  and  two-dimensional  modelling  would  appear  to  be  more  useful, 
as  discussed  in  Reference  14.  A  one-dimensional  model  may  be  useful,  but  only  in 
an  approximation  approach,  with  a  shaped  electric  field,  such  as  that  used  in 
Reference  21. 
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Figure  10.  Electric  Field  as  a  Function  of  y,  Distance  From  Anode,  Using 
Equations  ( 12)-(  16). 


34 


0.08 


4         6         8        10       12       14 
y/a 

Figure  1 1.  Species  Number  Density  as  a  Function  of  y ,  Distance  From  Anode,  Using 
Equations  (12)-(16). 
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E.      ANODE  FALL  VOLTAGE 

The  arc  discharge  operates  with  an  anode  voltage  drop  (Vaf)  which  should  be 
on  the  order  of  the  ionization  potential  of  the  gas  as  it  is  in  glow  discharges.  This 
voltage  drop  is  largely  non-ohmic  because  the  anode  region  is  typically  very  thin  and 
space-charge  controlled.  Since  a  significant  portion  of  the  anode  heating  comes  via 
the  oncoming  electrons  accelerating  through  the  fall  voltage,  it  is  of  interest  to 
minimize  this  voltage,  which  can  range  from  less  than  ten  to  over  40  volts  in  argon. 
The  question  arises  as  to  what  governs  this  anode  drop  and  why  is  it  such  a 
noticeable  function  of  current?  [Ref.  42] 

In  collisional  sheaths,  the  anode  voltage  drop  is  largely  governed  by  a  positive 
ion  generation  region  which  forms  in  front  of  the  anode.  The  production  of  ions 
reduces  the  space  charge  density  and  thereby  permits  operation  at  lower  voltages 
than  otherwise  possible.  Ions  are  most  often  created  by  electrons  within  their  last 
few  mean-free-paths  before  entering  the  anode;  these  ions  slowly  drift  away  from  the 
anode,  thereby  effectively  neutralizing  the  space  charge,  at  a  speed  proportional  to 
the  ratio  of  their  mobility  to  that  of  the  electrons.  At  the  anode,  the  electric  field 
is  a  maximum  and  the  electron  mean  energy  displays  a  corresponding  high. 

At  moderate  pressures,  the  sheath  is  very  thin  and  breakdown  can  be  visualized 
as  occurring  between  the  undisturbed  plasma  (which  acts  as  a  rich  source  of 
electrons)  and  the  anode  surface.  This  arrangement  has  the  attributes  that  represent 
"thermionic  arc  breakdown"  [Ref.  43],  a  source  of  electrons  which  is  independent  of 
the  breakdown  itself  and  relatively  small  spacings.   For  gases  that  allow  cumulative 
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ionization  (with  some  help  from  the  tail  of  the  Maxwellian  distribution  of  electron 
energies),  what  results  is  a  breakdown  voltage  appreciably  below  the  ionization 
potential.  This  then  could  be  an  explanation  for  the  low  voltage  breakdown 
observations  [Ref.  42].  Clearly,  gases  with  low  ionization  potentials  and  lots  of 
atomic  electron  energy  levels  are  preferred  (such  as  cesium  and  barium)  but  low- 
voltage  breakdown  has  been  observed  with  most  gases. 

The  increase  of  the  anode  fall  voltage  above  the  ionization  potential  has  been 
related  to  the  electron  Hall  parameter,  since  a  reduction  of  this  parameter  decreases 
that  voltage  and  corresponding  losses.  Control  of  the  local  magnetic  field  through 
the  use  of  and  array  of  permanent  magnets  as  well  as  implementation  of 
transpiration  cooling  (which  increases  the  electron  collision  frequency)  have  both 
yielded  some  encouraging  results.  Because  the  anode  fall  also  scales  up  with  J2/1^ 
it  is  conceivable  that  current  inhomogeneities  and  plasma  instabilities  which  are 
reflected  in  this  parameter  are  in  the  picture  as  well.  [Ref.  44] 

In  summary,  any  possible  reduction  of  the  anode  fall  voltage  will  hinge  on  a 
thorough  understanding  of  the  anode  region,  with  its  associated  sheath  and  ambipolar 
regions,  where  electron  temperature  effects,  ionization  effects,  and  magnetic  field 
effects  play  a  pivotal  role.  If  transpiration  cooling  is  present,  then  additional 
phenomena  of  fluid-dynamic  nature  may  come  into  play.  Experimental  observations 
with  atmospheric  discharges  indicate  the  possible  presence  of  convective  effects  at 
the  anode.  [Ref.  45] 
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IV.  TRANSPIRATION  COOLING 

Transpiration  cooling  of  the  anode  has  often  been  promoted  as  an  attractive 
means  of  recovering  a  large  portion  of  the  power  deposited  there.  Additionally,  the 
onset  of  melting  may  be  minimized  or  even  avoided  by  active  anode  cooling.  Rich, 
et.al.,  related  high  anode  temperatures  to  anode  spotting  [Ref.  37].  Similarly,  Park 
and  Choi  showed  that  low  thermal  diffusion  leads  to  erosion  and,  consequently, 
anode  damage  [Ref.  46].  Active  anode  cooling  via  transpiration  is  one  means  of 
ensuring  high  thermal  diffusion  and  extending  anode  lifetime.  Early  work  by 
Schoeck,  et.  al.,  [Ref.  47]  showed  that  up  to  80%  of  the  energy  deposited  in  the 
anode  is  recoverable  via  transpiration  cooling.  While  this  study  used  non-convective, 
high-intensity  argon  arcs,  it  is  reasonable  to  assume  that  this  effect  would  apply  in 
MPD  arcs  using  other  propellants.  Although  this  cooling  method  has  not  been 
studied  for  incorporation  in  plasma  thrusters  for  some  time,  it  has  been  recently 
considered  as  a  means  of  cooling  the  fuselage  of  the  National  Aerospace  Plane 
(NASP).  Plasma  thruster  designs  could  undoubtedly  gain  from  this  database,  and  due 
consideration  should  be  given  to  this  cooling  approach  for  the  anode. 

For  a  given  mass  transfer  flow  rate,  the  heat  flux  reduction  to  a  surface  is 
inversely  proportional  to  the  molecular  weight  of  the  injected  gas.  Use  of  the 
propellant  as  coolant  as  well  as  fuel  would  eliminate  the  need  for  additional  tankage 
and  pumps,  simplifying  the  design  considerably.   Lithium  has  been  considered  to  be 


38 


the  propellant  of  choice,  primarily  because  of  its  low  molecular  weight,  its  favorable 

ionization  potential,  and  its  low-volume  tankage  properties.  It  has  a  relatively  low 

first  ionization  potential  of  5.4  eV  and  a  high  second  ionization  potential  of  75.6  eV. 

This  single  ionization  potential  range  of  over  70  eV  compares  to  approximately  20 

eV  for  Cesium  and  27  e V  for  Potassium  [Ref.  48].  This  provides  a  broad  temperature 

range  within  which  only  single  ionization  will  occur.  Large  temperatures  must  be 

reached  within  the  gas  before  double  ionization  occurs.  As  the  gas  temperature  is 

increased  several  thousand  degrees  Kelvin,  it  undergoes  ionization  and  disassociation. 

Thermal  energy  deposited  can  be  recovered  through  nozzle  expansion  at  the  exit. 

However,  residence  time  of  the  gas  is  not  long  enough  to  ensure  recombination. 

Thus,  the  energy  invested  in  ionization  and  disassociation  will  be  lost.  [Ref.  49] 

Lithium  has  been  shown  to  produce  specific  impulse  figures  in  excess  of  7,000 

seconds  at  70%  efficiency  in  steady-state  thrusters,  [Ref.  50]  whereas  all  other 

propellants  have  been  limited  to  less  than  3,000  seconds  specific  impulse  at  less  than 

40%  efficiency  [Ref.  6].   Subramaniam  has  concluded  that: 

...regenerative  cooling  of  anodes  (at  the  specific  impulse  values  in  the  MPD 
regime)  is  possible  only  with  hydrogen  or  with  alkali-metal  propellants,  notably 
lithium.  In  the  latter  case,  the  ideal  anode  operating  mode  would  be 
evaporation  and  ionization  of  the  propellant  on  the  porous  or  wetted  anode 
surface,  resulting  in  increased  ion  current  fraction,  reduced  anode  voltage  fall 
and  utilization  of  part  of  the  anode  loss  energy  [Ref.  51]. 

Liquid  coolants,  as  well  as  reducing  storage  requirements,  offer  the  advantage  of 

providing  latent  heat  of  vaporization  for  energy  disposal.  However,  design  problems 

can  occur  if  the  liquid  is  allowed  to  vaporize  within  the  porous  structure.   Problems 
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arise  due  to  the  abrupt  increase  in  pressure  gradient  as  the  coolant  vaporizes.  Since 
coolant  flow  generally  have  three-dimensional  characteristics,  the  flow  will  be 
diverted  around  the  vapor  bubbles  and  hot  spots  often  develop.  The  technical 
practicality  of  using  molten  lithium  to  cool  a  porous  tungsten  anode  would  seem  to 
be  beyond  current  technology.  On  the  other  hand,  the  products  of  decomposition  of 
hydrazine  (gaseous  hydrogen  and  ammonia)  have  proven  to  be  efficient  and  practical 
coolants  [Ref.  52]. 

Given  the  performance  figures  above,  using  an  auxiliary  coolant  gas  even  with 
high  molecular  weight  (e.g.,  NH3,  N2,  CH4,  etc.)  which  could  serve  as  a  propellant 
once  released  from  a  porous  hot  tungsten  anode  surface  would  seem  to  more 
practical,  vice  dealing  with  molten  lithium.  Experimental  studies  would  be  needed 
to  compare  the  approaches.  Kuriki  and  Suzuki  performed  experiments  with  a  quasi- 
steady  MPD  thruster  to  study  the  effect  of  anode  gas  injection  (Argon).  At  high 
currents  of  up  to  10  kA,  increases  in  thrust,  specific  impulse,  and  flow  discharge 
stability  were  observed  [Ref.  53]. 

There  is  some  question  as  to  the  likelihood  of  current  constriction  resulting 
from  anode  gas  injection.  In  such  a  case,  swirling  or  circulating  the  propellant  gas 
would  help  to  move  any  footpoints  that  developed  around  the  anode  surface  and 
prevent  them  from  becoming  fully-developed  spots.  Additionally,  an  applied 
magnetic  field  could  serve  to  circulate  the  footpoints  as  well.  The  unique  advantage 
of  transpiration  cooling  hinges  on  providing  effective  anode  cooling  while  supplying 
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hot  propellant,  but  the  real  benefit  will  depend  on  how  small  the  amount  of  coolant 
required  really  will  be. 

Transpiration  cooling  has  proven  to  be  as  desirable  as  it  is  challenging.  It  is 
complicated  to  implement,  with  associated  reliability  problems  and  difficulty  of 
analytical  predictions.  While  the  production  of  thicker  boundary  layers  is  largely 
ineffective  against  the  electron  flux  heating,  the  cooling  itself  is  most  efficient  and  a 
substantial  fraction  of  the  energy  transferred  to  the  anode  is  recoverable.  The 
arguments  of  Chapter  Three  indicate  that  a  reduction  of  the  electron  temperature 
in  the  anode  would  have  the  desirable  effect  of  reducing  the  initial  current  spotting 
which  can  be  conjectured  to  be  the  path  that  leads  to  anode  arc  spots.  This  electron 
temperature  reduction  can  be  done  most  effectively  by  polyatomic  gases  (which  have 
a  high  5-loss  factor)  emanating  from  the  anode  surface  [Ref.  54]. 

The  arguments  relating  to  transpiration  cooling  might  be  summarized  as 
follows: 
Favorable  Outcomes 

•  No  separate  cooling  mechanism  for  anode  required, 

•  Adds  "hot"  propellant  to  exhaust  "recovering"  most  of  the  electrical  power  loss  to 
the  anode, 

•  Quenches  Te  thus  likely  to  postpone  anode  spotting  and  reduce  the  heating 
associated  with  the  electron  thermal  energy  (5kTe/2e), 

•  Reduces  bulk  convective  heating, 

•  Reduces  the  local  electron  Hall  parameter  by  increasing  the  collision  frequency, 
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Favorable  Outcomes  (cont'd.') 

•  Allows  for  some  radiation  cooling  from  the  hot  tungsten  surface  (about  120 
watts/cm2  at  2800°K  [Ref.  51]), 

Hydrogen/ammonia  gases   flowing   through   hot  porous   (sintered)   tungsten 
represent  a  compatible,  proven  technology. 

Unfavorable  Outcomes 

•  May  decrease  the  electrical  conductivity  in  the  anode  region, 

•  May  destabilize  the  ionization  processes  in  the  sheath  and  bring  about  significant 
fluctuations  in  the  current, 

•  Disrupts  "cathode  jet"  in  front  of  the  anode  with  unpredictable  consequences, 

•  Introduces  propellant  which  may  not  be  hot  enough,  not  ionized  enough,  or  not 
in  the  proper  place  for    JxB  acceleration, 

•  Transpiration  cooling  through  a  porous  (tungsten)  anode  is  a  difficult  design 
problem. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 

Plasma  thrusters  offer  distinct  advantages  in  terms  of  payload  delivered  for 
interplanetary  missions,  as  well  as  for  orbital  transfer.  A  recent  comparison 
completed  by  Choueiri,  Kelly,  and  Jahn  shows  a  mass  savings  of  65  tons  for  an 
orbital  transfer  from  low  Earth  orbit  to  geosynchronous  Earth  orbit  using  a 
quasisteady  MPD  thruster  as  opposed  to  an  advanced  chemical  thruster.6  This 
superior  performance  comes  at  the  expense  of  low  thrust-to-weight  ratio  and  long 
transit  time.  However,  given  the  large  cargo/logistic  requirements  of  a  manned 
interplanetary  mission,  delivery  of  payload  must  be  maximized.  Thus,  further  work 
to  characterize  more  fully  thruster  behavior  and  anode  contributions  in  particular  are 
certainly  warranted.  [Ref.  55] 

The  "cut-and-paste"  method  used  to  generate  Figures  10  and  11  is  not  of 
practical  use  as  a  modelling  method,  due  to  the  large  effort  involved.  It  did  produce 
the  expected  electric  field  and  species  number  and  current  density  plots  near  the 
anode,  but  failed  to  produce  the  entire  sheath  out  to  the  ambipolar  region.  The 
nonlinearity  of  the  equation  set  led  to  a  quickly  deteriorating  solution.  A  more 
practical  approach  using  nonlinear  dynamics  and/or  chaos  must  be  developed  to 
model  the  sheath  numerically. 


6This  assumes  a  specific  impulse  of  2,000  seconds,  600  kW  of  input  power,  and 
a  270-day  transit  time. 
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Depending  on  the  propellant  mass  fraction  used  for  cooling,  the  transpiration 
scheme  discussed  above  presents  some  rather  unique  advantages.  A  hot  anode 
which  uses  only  a  small  amount  of  propellant  for  cooling  need  not  be  penalized  for 
any  lost  thrust.  If  in  addition,  we  increase  anode  lifetime  by  delaying  the  formation 
of  anode  arc  spots,  then  the  scheme  is  all  the  more  desirable.  A  decrease  of  the 
electron  temperature  in  the  vicinity  of  the  anode  may  bring  about  a  more 
homogeneous  flow  of  current  and  a  reduction  in  the  heating  effect  associated  with 
the  high  electron  kinetic  energy.  Recovery  of  the  heat  deposited  at  the  anode  would 
be  most  important  if  the  propellant  fraction  is  high.  In  such  case,  nozzle  expansion 
of  the  hot-propellant/coolant-gas  might  be  implemented. 

Means  of  limiting  anode  losses  through  decreasing  anode  fall  voltages  were 
discussed,  including  the  control  of  the  local  Hall  parameter  and  the  implementation 
of  thermionic  arc  breakdown.  The  electrical  conductivity  (of  a  nonreacting  plasma) 
could  possibly  decrease  as  a  result  of  transpiration  cooling  and  this  might  increase 
the  anode  fall  voltage. 

Additional  work  needs  to  be  done  in  the  following  areas: 

•  Investigate  effectiveness  of  nonlinear  dynamics  and  chaos  in  solving  sheath 
equation  set, 

•  Incorporate  adequate  one-  or  two-dimensional  sheath  modelling  in  quasisteady 
MPD  numerical  codes, 
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•  Investigate  the  role  that  fluid  dynamic  effects  play  in  MPD  thruster  anode 
discharges, 

•  Investigate  the  effect  of  transpiration  cooling  on  current  and  plasma  stability,  as 
well  as  on  thruster  performance  and  lifetime, 

•  Determine  effectiveness  of  transpiration  cooling's  increase   of  the  collision 
frequency  parameter, 

•  Compare  performance  of  gaseous  propellant/coolants  versus  hybrid  designs  with 
lithium  propellant/gaseous  coolant, 

•  Determine  if  required  percentage  of  propellant  gas  as  coolant  is  practical  (e.g.,  less 
than  10%), 

•  Investigate  effect  of  surface  imperfections  as  focal  points  for  current  constrictions 
and  as  precursors  to  anode  spotting. 
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APPENDIX  A 


The  following  software  includes  the  calling  program,  SHEATH,  its  two 
subroutines,  FCNJ  and  YDOT,  and  the  DGEAR  integrator.  The  latter  is  quite 
extensive  in  length  and  includes  ten  subroutines,  including  the  following:  DGRST, 
DGRCS,  DGRPS,  DGRIN,  LUDATF,  LUELMF,  LEQT1B,  UERTST,  UGETIO,  and  USPKD .  A 
detailed  discussion  may  be  found  in  IMSL  literature  or  Reference  39 . 

it******************************************************* 

Program  Sheath  000010 

C  000020 

C Calling  program  for  DGEAR  integrator.  Initial  conditions  are  000030 

C     input  via  READ  statements  and  keyboard  entry.  Output  is  to  data    000040 

C      files  via  the  DGRST  subroutine.  Diagnostic  check  of  output  via  000050 

C      Figure  4  printed  to  data  file  from  DGRST  subroutine.  Consult  000060 

C      DGEAR  comments  for  variable  descriptions  not  listed  below.  000061 

C 000062 

REAL  E,K,EPS,TI,EF,EFINF,DI,DE,NUINF,C1,K1,A  000070 

INTEGER  N, METH, MITER, INDEX, IWK(l) , IER, STEP  000080 

REAL*8  X,H, Y(5) , XEND , TOL , WK  000090 

EXTERNAL  YDOT,  FCNJ,  DGEAR  000100 

COMMON/ CONST /E , K, EPS , TI , EF , EFINF , NINF , DI , DE , VINF ,C1,K1,A  000110 

C  000120 

C Constants  000121 

C      CI  and  Kl  are  constants  describing  the  ionization  coefficient.  000122 

C     They  are  taken  from  the  data  plotted  in  Figure  8.  The  000123 
C     coefficient  is  equal  to  the  nondimensionalized  electric  potential  000124 

C     raised  to  the  Kl  power  and  then  multiplied  by  CI.  000125 

C     In  this  way,  the  ionization  coefficient  is  allowed  to  vary  in  000126 

C     proportion  to  the  strength  of  the  electric  potential.  000127 

C 000128 

WRITE (*,*)' Input  value  for  CI  (format  6E3):'  000129 

READ(*,*)C1  000130 

WRITE (*,*)C1  000131 

WRITE (*,*)' Input  value  for  Kl  (format  6E3):'  000132 

READ(*,*)K1  000133 

WRITE(*,*)K1  000134 
C      Initial  conditions  for  species  number  density,  electric  potential  000135 

C     and  species  current  density  are  now  input  (ni , ne, E, je, j i) .  000136 

C--- 000137 

WRITE (*,*)' Input  values  for  y(l)  through  y(5)  (format  5(6E3)):'    000138 

READ(*,*)y(l)  ,y(2),y(3),y(4),y(5)  000139 

WRITE  (*,*)  y(l)  ,y(2),y(3),y(4),y(5)  00014  0 

C      Following  constants  are  for  plasma  described  in  Reference  21  000141 

C      (6,000  K,  Init  E=20,000  V/m,  Final  E=l,200  V/m)  000142 

E=1.6E-19  000143 

K=1.38E-23  000150 

EPS=8 .854E-12  000160 

TI=6E3  000170 

EF=2E5  000180 

EFINF=1.2E4  000190 

DI=1.724E-4  000200 

DE=1.724E-1  000210 

VIO  =  2 .E6  000215 

VINF  =  4 .93E-7  000220 

C  000230 
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C      A  is  plasma  characteristic  length  which  shows  potential  drop. 

000240 

C      A  =  (  (EPS*EF) / (E*NINF) )  =  1.107E-6 

A  =  1.107E-5  000250 

X  =  0.01  000270 

XEND  =  10.  000280 

H  =  le-6  000290 

TOL  =  1E-6  000300 

METH  =  2  000305 

MITER  =  1  000310 

INDEX  =  1  000320 

N=5  000330 

IWK(l)  =  5  000340 

WK  =  18000.  000350 

IER  =  0  000360 

OPEN (UNIT=8,FILE=' SHEATH. DAT' , STATUS =' UNKNOWN' )  000370 
CALL  DGEAR2 (N, YDOT, FCNJ, X, H, Y, XEND , TOL, METH, MITER, INDEX, IWK,WK,    000380 

+IER,STEP)  000390 

DO  3  1=0, N  000400 

DO  2  J=0,100  000410 

WRITE(*,*) J,Y(I)  000420 

WRITE (8,1) J,Y(I)  000430 

1  FORMAT (T2,F5. 1,5 (5X,D9.2) )  000440 

2  CONTINUE  000450 

3  CONTINUE  000460 
WRITE (*,*) 'Total  Steps  =  ', STEP, 'Final  Step  Size  =  ',H,            000470 

+' Error  Code  =  ' , IER  000480 

CLOSE (UNIT=8)  000490 

STOP  000500 

END  000510 

C  DUMMY  SUBROUTINE  FCNJ 
p*********************************** 

SUBROUTINE  FCNJ (N, X, Y, PD)  1 

INTEGER  N  2 

REAL  Y(N)  ,PD(N,N)  ,X  .         3 

RETURN  4 

END  5 

C  SUBROUTINE  YDOT 
p*********************************** 

SUBROUTINE  YDOT (N, X , Y , YPRIME , eprime , eprime2 ) 

REAL* 8  X ,  Y ( 5 )  , YPRIME ( 5 )  , NUI , eprime , epr ime2 
REAL  E , K, EPS , TI , EF , EFINF , NINF , DI , DE , VINF , CI , Kl , A, Bl , B2 , B3 , B4 
COMMON/CONST/E , K , EPS , TI , EF , EFINF , NINF , DI , DE , VIO , VINF , CI , Kl , A 

VI  =  CI  *  (Y(3) **K1) 

VIT  =  VI  /  VIO 
C     Following  constants  are  the  bracketed  values  in  Equations  12-16. 
C     A  is  left  as  a  variable. 

C - 

C      Bl  =  ( (E*EPS)/(K*TI) )  *  A 

Bl=  3 .86E5  *  A 
C      B2  =  ( (E*EFINF) /(K*TI) )  *  A 

B2  =  2.32E4  *  A 
C      B3  =  ( (E*NINF) / (EF*EPS) )  *  A 

B3  =  9 .04E5  *  A 
C      B4  =  ( (VINF*K*TI) / (E*DE*EFINF) )  *  A 

B4  =  2  .62E-21  *  A 
C     Alpha  =  2 -body  recombination  coefficient  (fm.  Laser  Kinetics 
C     Handbook  (AFWL-TR-74 -216 ,  1974))  (cm3/sec) 

Alpha  =  9 . e-8 
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c-  - 

C  FIVE  FIRST  ORDER  EQUATIONS  -  Equations  12-16 

C 

C     Ni 

YPRIME(l)  =   (B  *  Y(l)  *  Y(3))  -  Y(5) 
C      Ne 

YPRIME(2)  =  -  (B  *  Y(2)  *  Y(3))  +  Y(4) 
C      E 

YPRIMEO)  =   B3  *  (Y(l)  -  Y(2)) 
C      je 

YPRIMEU)  =  -B4  *  Y(2)  *  (VIT  -  (ALPHA  *  Y(l))) 
C      ji 

YPRIME(5)  =   B4  *  Y(2)  *  (VIT  -  (ALPHA  *  Y(l))) 
C 

C  —  Diagnostic  Check  of  first, second  derivatives 

C 

eprime  =  y(l)  -  y(2) 

eprime2  =  yprime(l)  -  yprime(2) 
C 

RETURN 

END 
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IMSL  ROUTINE  NAME 
•modified  to  return  # 
COMPUTER 
LATEST  REVISION 
PURPOSE 

USAGE 

ARGUMENTS     N 

FCN 


-  DGEAR  DGEA0010 

DGEA002  0 
of  steps  via  variable  "step"  in  subroutine  call  + 

DGEA004  0 

-  IBM/DOUBLE  DGEA0050 

DGEA0060 

-  NOVEMBER  1,  1984  DGEA0070 

DGEA0080 
DIFFERENTIAL  EQUATION  SOLVER  -  VARIABLE  ORDER  DGEA0  09  0 


FCNJ 


ADAMS  PREDICTOR  CORRECTOR  METHOD  OR 
GEARS  METHOD 

CALL  DGEAR  (N, FCN, FCNJ, X, H, Y,XEND, TOL,  METH, 
MITER, INDEX, IWK, WK, IER) 

INPUT  NUMBER  OF  FIRST-ORDER  DIFFERENTIAL 

EQUATIONS . 
NAME  OF  SUBROUTINE  FOR  EVALUATING  FUNCTIONS. 
(INPUT) 


DGEA0100 
DGEA0110 
DGEA0120 
DGEA013  0 
DGEA014  0 
DGEA0150 
DGEA0160 
DGEA0170 
DGEA0180 
DGEA019  0 


THE  SUBROUTINE  ITSELF  MUST  ALSO  BE  PROVIDED  DGEA0200 


BY  THE  USER  AND  IT  SHOULD  BE  OF  THE 
FOLLOWING  FORM 

SUBROUTINE  FCN  (N,X, Y, YPRIME) 

REAL  X,Y(N) , YPRIME (N) 


DGEA0210 

DGEA0220 

DGEA023  0 

DGEA024  0 

DGEA0250 

DGEA0260 

DGEA0270 

FCN  SHOULD  EVALUATE  YPRIME (1) , . . . , YPRIME (N)  DGEA0280 

GIVEN  N,X,  AND  Y(l) , . . . ,Y(N) .  YPRIME (I)    DGEA0290 

IS  THE  FIRST  DERIVATIVE  OF  Y(I)  WITH       DGEA0300 

RESPECT  TO  X.  DGEA0310 

FCN  MUST  APPEAR  IN  AN  EXTERNAL  STATEMENT  IN  DGEA0320 

THE  CALLING  PROGRAM  AND  N,X, Y (1) , . . . , Y (N)  DGEA0330 

MUST  NOT  BE  ALTERED  BY  FCN.  DGEA034  0 

NAME  OF  THE  SUBROUTINE  FOR  COMPUTING  THE       DGEA0350 

JACOBIAN  MATRIX  OF  PARTIAL  DERIVATIVES.      DGEA03  60 

(INPUT)  DGEA037  0 

THE  SUBROUTINE  ITSELF  MUST  ALSO  BE  PROVIDED  DGEA03  80 


BY  THE  USER. 
IF  MITER=1  IT  SHOULD  BE  OF  THE  FOLLOWING 
FORM 

SUBROUTINE  FCNJ  (N,X,Y,PD) 
REAL  X,Y(N)  ,PD(N,N) 


FCNJ  MUST  EVALUATE  PD(I,J) ,  THE  PARTIAL 
DERIVATIVE  OF  YPRIME (I)  WITH  RESPECT  TO 
Y(J) ,  FOR  1=1, N  AND  J=1,N. 
IF  MITER=  -1  IT  SHOULD  BE  OF  THE  FOLLOWING 
FORM 

SUBROUTINE  FCNJ  (N,X,Y,PD) 
REAL  X,Y(N)  ,PD(1) 


DGEA039  0 

DGEA04  00 

DGEA0410 

DGEA04  20 

DGEA0430 

DGEA044  0 

DGEA04  50 

DGEA04  60 

DGEA04  70 

DGEA04  80 

DGEA0490 

DGEA0500 

DGEA0510 

DGEA052  0 

DGEA0530 

DGEA054  0 

FCNJ  MUST  EVALUATE  PD  IN  BAND  STORAGE  MODE.  DGEA0550 

THAT  IS,  PD(N*  (J-I+NLO+I)  IS  THE  PARTIAL  DGEA0560 

DERIVATIVE  OF  YPRIME (I)  WITH  RESPECT  TO    DGEA0570 

Y(J)  .   NLC  IS  THE  NUMBER  OF  LOWER  DGEA0580 

CODIAGONALS  FOR  THE  BAND  MATRIX.  DGEA0590 

FCNJ  MUST  APPEAR  IN  AN  EXTERNAL  STATEMENT  INDGEA0600 

THE  CALLING  PROGRAM  AND  N, X, Y (1) , . . . , Y (N)  DGEA0610 
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H 


XEND 


TOL 


METH 


MITER 


MUST  NOT  BE  ALTERED  BY  FCNJ .  DGEA0620 

FCNJ  IS  USED  ONLY  IF  MITER  IS  EQUAL  TO  DGEA063  0 

1  OR  -1.  OTHERWISE  A  DUMMY  ROUTINE  CAN  DGEA0640 

BE  SUBSTITUTED.  SEE  REMARK  1.  DGEA0650 

INDEPENDENT  VARIABLE.  (INPUT  AND  OUTPUT)  DGEA0660 

ON  INPUT,  X  SUPPLIES  THE  INITIAL  VALUE  DGEA067  0 

AND  IS  USED  ONLY  ON  THE  FIRST  CALL.  DGEA0680 

ON  OUTPUT,  X  IS  REPLACED  WITH  THE  CURRENT  DGEA0690 

VALUE  OF  THE  INDEPENDENT  VARIABLE  AT  WHICHDGEA0700 

INTEGRATION  HAS  BEEN  COMPLETED.  DGEA0710 

INPUT/OUTPUT.  DGEA0720 

ON  INPUT,  H  CONTAINS  THE  NEXT  STEP  SIZE  IN  DGEA0730 

X.  H  IS  USED  ONLY  ON  THE  FIRST  CALL.  DGEA074  0 

ON  OUTPUT,  H  CONTAINS  THE  STEP  SIZE  USED  DGEA0750 

LAST,  WHETHER  SUCCESSFULLY  OR  NOT.  DGEA0760 

DEPENDENT  VARIABLES,  VECTOR  OF  LENGTH  N.  DGEA0770 

(INPUT  AND  OUTPUT)  DGEA0780 

ON  INPUT,  Y(l)  ,  .  .  . ,Y(N)  SUPPLY  INITIAL  DGEA0790 

VALUES.  DGEA0800 

ON  OUTPUT,  Y(l) , . . . ,Y(N)  ARE  REPLACED  WITH  DGEA0810 

A  COMPUTED  VALUE  AT  XEND.  DGEA082  0 

INPUT  VALUE  OF  X  AT  WHICH  SOLUTION  IS  DESIRED  DGEA083  0 

NEXT.  INTEGRATION  WILL  NORMALLY  GO  DGEA084  0 

BEYOND  XEND  AND  THE  ROUTINE  WILL  INTERPOLATEDGEA0850 

TO  X  =  XEND.  DGEA0860 

NOTE  THAT  (X-XEND)  *H  MUST  BE  LESS  THAN  DGEA0870 

ZERO  (X  AND  H  AS  SPECIFIED  ON  INPUT) .  DGEA0880 

INPUT  RELATIVE  ERROR  BOUND.  TOL  MUST  BE  DGEA089  0 

GREATER  THAN  ZERO.  TOL  IS  USED  ONLY  ON  THE  DGEA0900 

FIRST  CALL  UNLESS  INDEX  IS  EQUAL  TO  -1.  DGEA0910 

TOL  SHOULD  BE  AT  LEAST  AN  ORDER  OF  DGEA0920 

MAGNITUDE  LARGER  THAN  THE  UNIT  ROUNDOFF  DGEA0930 

BUT  GENERALLY  NOT  LARGER  THAN  .001.  DGEA094  0 

SINGLE  STEP  ERROR  ESTIMATES  DIVIDED  BY  DGEA0950 

YMAX(I)  WILL  BE  KEPT  LESS  THAN  TOL  IN  DGEA0960 

ROOT -MEAN -SQUARE  NORM  (EUCLIDEAN  NORM  DGEA0970 

DIVIDED  BY  SQRT(N))  .  THE  VECTOR  YMAX  OF  DGEA0980 

WEIGHTS  IS  COMPUTED  INTERNALLY  AND  STORED  DGEA0990 

IN  WORK  VECTOR  WK.  INITIALLY  YMAX  (I)  IS  DGEA1000 

THE  ABSOLUTE  VALUE  OF  Y(I)  ,  WITH  A  DEFAULT  DGEA1010 

VALUE  OF  ONE  IF  Y(I)  IS  EQUAL  TO  ZERO.  DGEA1020 

THEREAFTER,  YMAX (I)  IS  THE  LARGEST  VALUE  DGEA1030 

OF  THE  ABSOLUTE  VALUE  OF  Y(I)  SEEN  SO  FAR,  DGEA1040 

OR  THE  INITIAL  VALUE  OF  YMAX  (I)  IF  THAT  IS  DGEA1050 

LARGER.  DGEA1060 

INPUT  BASIC  METHOD  INDICATOR.  DGEA1070 

USED  ONLY  ON  THE  FIRST  CALL  UNLESS  INDEX  IS  DGEA1080 

EQUAL  TO  -1.  DGEA109  0 

METH  =  1,  IMPLIES  THAT  THE  ADAMS  METHOD  IS  DGEA1100 

TO  BE  USED.  DGEA1110 

METH  =  2,  IMPLIES  THAT  THE  STIFF  METHODS  OF  DGEA1120 

GEAR,  OR  THE  BACKWARD  DIFFERENTIATION  DGEA1130 

FORMULAE  ARE  TO  BE  USED.  DGEA1140 

INPUT  ITERATION  METHOD  INDICATOR.  DGEA1150 

MITER  =  0,  IMPLIES  THAT  FUNCTIONAL  DGEA1160 

ITERATION  IS  USED.  NO  PARTIAL  DGEA1170 

DERIVATIVES  ARE  NEEDED.  A  DUMMY  FCNJ  DGEA1180 

CAN  BE  USED.  DGEA119  0 

MITER  =  1,  IMPLIES  THAT  THE  CHORD  METHOD  DGEA1200 

IS  USED  WITH  AN  ANALYTIC  JACOBIAN.  FOR  DGEA1210 

THIS  METHOD,  THE  USER  SUPPLIES  DGEA1220 

SUBROUTINE  FCNJ.  DGEA1230 
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INDEX 


IWK 
WK 


MITER  =  2,  IMPLIES  THAT  THE  CHORD  METHOD  DGEA124  0 

IS  USED  WITH  THE  JACOB  IAN  CALCULATED  DGEA1250 

INTERNALLY  BY  FINITE  DIFFERENCES.  DGEA1260 

A  DUMMY  FCNJ  CAN  BE  USED.  DGEA1270 

MITER  =  3,  IMPLIES  THAT  THE  CHORD  METHOD  DGEA1280 

IS  USED  WITH  THE  JACOB  IAN  REPLACED  BY  DGEA1290 

A  DIAGONAL  APPROXIMATION  BASED  ON  A  DGEA13  00 

DIRECTIONAL  DERIVATIVE.  DGEA1310 

A  DUMMY  FCNJ  CAN  BE  USED.  DGEA132  0 

MITER  =  -1  OR  -2,  IMPLIES  USE  THE  SAME  DGEA133  0 

METHOD  AS  FOR  MITER=  1  OR  2,  RESPECTIVELY, DGEA1340 

BUT  USING  A  BANDED  JACOBIAN  MATRIX.   IN  DGEA1350 

THESE  TWO  CASES  BANDWIDTH  INFORMATION  DGEA13  60 

MUST  BE  PASSED  TO  DGEAR  THROUGH  THE  DGEA13  70 

COMMON  BLOCK  DGEA13  80 

COMMON  /DBAND/  NLC,NUC  DGEA1390 

WHERE  NLC=NUMBER  OF  LOWER  CODIAGONALS  DGEA14  00 

NUC=NUMBER  OF  UPPER  CODIAGONALS  DGEA1410 

INPUT  AND  OUTPUT  PARAMETER  USED  TO  INDICATE  DGEA14  2  0 

THE  TYPE  OF  CALL  TO  THE  SUBROUTINE.   ON  DGEA14  3  0 

OUTPUT  INDEX  IS  RESET  TO  0  IF  INTEGRATION  DGEA144  0 

WAS  SUCCESSFUL.   OTHERWISE,  THE  VALUE  OF  DGEA1450 

INDEX  IS  UNCHANGED.  DGEA14  60 

ON  INPUT,  INDEX  =  1,  IMPLIES  THAT  THIS  IS  THE  DGEA1470 


DGEA14  80 

IS  NOT  DGEA1490 

DGEA1500 


FIRST  CALL  FOR  THIS  PROBLEM. 
ON  INPUT,  INDEX  =  0,  IMPLIES  THAT  THIS 

THE  FIRST  CALL  FOR  THIS  PROBLEM. 

ON  INPUT,  INDEX  =  -1,  IMPLIES  THAT  THIS  IS  NOTDGEA1510 

THE  FIRST  CALL  FOR  THIS  PROBLEM,  AND  THE     DGEA1520 

USER  HAS  RESET  TOL .  DGEA1530 

ON  INPUT,  INDEX  =  2,  IMPLIES  THAT  THIS  IS  NOT  DGEA1540 

THE  FIRST  CALL  FOR  THIS  PROBLEM.  INTEGRATIONDGEA1550 

IS  TO  CONTINUE  AND  XEND  IS  TO  BE  HIT  EXACTLYDGEA1 5 6 0 

(NO  INTERPOLATION  IS  DONE) .  THIS  VALUE  OF    DGEA1570 

INDEX  ASSUMES  THAT  XEND  IS  BEYOND  THE        DGEA1580 

CURRENT  VALUE  OF  X.  DGEA159  0 

ON  INPUT,  INDEX  =  3,  IMPLIES  THAT  THIS  IS  NOT  DGEA1600 

THE  FIRST  CALL  FOR  THIS  PROBLEM.  INTEGRATIONDGEA1610 

IS  TO  CONTINUE  AND  CONTROL  IS  TO  BE  RETURNEDDGEA162  0 

TO  THE  CALLING  PROGRAM  AFTER  ONE  STEP.  XEND  DGEA1630 

IS  IGNORED.  DGEA164  0 

INTEGER  WORK  VECTOR  OF  LENGTH  N.  USED  ONLY  IF  DGEA1650 

MITER  =  1  OR  2  DGEA1660 

REAL  WORK  VECTOR  OF  LENGTH  4*N+NMETH+NMITER. 

THE  VALUE  OF  NMETH  DEPENDS  ON  THE  VALUE  OF 

METH. 

IF  METH  IS  EQUAL  TO  1, 
NMETH  IS  EQUAL  TO  N*13 . 


IF  METH  IS  EQUAL  TO  2, 
NMETH  IS  EQUAL  TO  N*6. 


DGEA167  0 
DGEA1680 
DGEA1690 
DGEA1700 
DGEA1710 
DGEA1720 
DGEA173  0 


THE  VALUE  OF  NMITER  DEPENDS  ON  THE  VALUE  OF  DGEA174  0 

MITER.  DGEA175  0 

IF  MITER  IS  EQUAL  TO  1  OR  2 ,  DGEA1760 

NMITER  IS  EQUAL  TO  N*(N+1)  DGEA1770 

IF  MITER  IS  EQUAL  TO  -1  OR  -2,  DGEA1780 

NMITER  IS  EQUAL  TO  (2*NLC+NUC+3 ) *N  DGEA1790 

WHERE  NLC=NUMBER  OF  LOWER  CODIAGONALS  DGEA1800 

NUC=NUMBER  OF  UPPER  CODIAGONALS  DGEA1810 

IF  MITER  IS  EQUAL  TO  3,  DGEA1820 

NMITER  IS  EQUAL  TO  N.  DGEA1830 

IF  MITER  IS  EQUAL  TO  0,  DGEA184  0 

NMITER  IS  EQUAL  TO  1.  DGEA1850 
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IER 


PRECI 


Step 
S ION/HARDWARE   - 


WK  MUST  REMAIN  UNCHANGED  BETWEEN  SUCCESSIVE  DGEA1860 

CALLS  DURING  INTEGRATION.  DGEA1870 

ERROR  PARAMETER.  (OUTPUT)  DGEA1880 

WARNING  ERROR  DGEA1890 

IER  =  33,  IMPLIES  THAT  X+H  WILL  EQUAL  X  ON  JGEA1900 

THE  NEXT  STEP.  THIS  CONDITION  DOES  NOT  DGEA1910 

FORCE  THE  ROUTINE  TO  HALT.  HOWEVER,  IT  DGEA1920 

DOES  INDICATE  ONE  OF  TWO  CONDITIONS.  DGEA1930 

THE  USER  MIGHT  BE  REQUIRING  TOO  MUCH  DGEA194  0 

ACCURACY  VIA  THE  INPUT  PARAMETER  TOL.  DGEA1950 

IN  THIS  CASE  THE  USER  SHOULD  CONSIDER  DGEA1960 

INCREASING  THE  VALUE  OF  TOL.  THE  OTHER  DGEA1970 

CONDITION  WHICH  MIGHT  GIVE  RISE  TO  THIS  DGEA1980 

ERROR  MESSAGE  IS  THAT  THE  SYSTEM  OF  DGEA1990 

DIFFERENTIAL  EQUATIONS  BEING  SOLVED  DGEA2000 

IS  STIFF  (EITHER  IN  GENERAL  OR  OVER  DGEA2010 

THE  SUB INTERVAL  OF  THE  PROBLEM  BEING  DGEA2  020 

SOLVED  AT  THE  TIME  OF  THE  ERROR) .  IN  DGEA203  0 

THIS  CASE  THE  USER  SHOULD  CONSIDER  DGEA2  04  0 

USING  A  NONZERO  VALUE  FOR  THE  INPUT  DGEA2050 

PARAMETER  MITER.  DGEA2060 

WARNING  WITH  FIX  ERROR  DGEA2  070 

IER  =  66,  IMPLIES  THAT  THE  ERROR  TEST  DGEA2  080 

FAILED.  H  WAS  REDUCED  BY  .1  ONE  OR  MORE  DGEA2090 

TIMES  AND  THE  STEP  WAS  TRIED  AGAIN  DGEA2100 

SUCCESSFULLY.  DGEA2110 

IER  =67,  IMPLIES  THAT  CORRECTOR  DGEA2120 

CONVERGENCE  COULD  NOT  BE  ACHIEVED.  DGEA2130 

H  WAS  REDUCED  BY  .1  ONE  OR  MORE  TIMES  AND  DGEA214  0 

THE  STEP  WAS  TRIED  AGAIN  SUCCESSFULLY.  DGEA2150 

TERMINAL  ERROR  DGEA2160 

IER  =  132,  IMPLIES  THE  INTEGRATION  WAS  DGEA2170 

HALTED  AFTER  FAILING  TO  PASS  THE  ERROR  DGEA2180 

TEST  EVEN  AFTER  REDUCING  H  BY  A  FACTOR  DGEA2190 

OF  1.0E10  FROM  ITS  INITIAL  VALUE.  DGEA2200 

SEE  REMARKS.  DGEA2210 

IER  =  133,  IMPLIES  THE  INTEGRATION  WAS  DGEA2220 

HALTED  AFTER  FAILING  TO  ACHIEVE  DGEA2230 

CORRECTOR  CONVERGENCE  EVEN  AFTER  DGEA224  0 

REDUCING  H  BY  A  FACTOR  OF  1.0E10  FROM  DGEA2250 

ITS  INITIAL  VALUE.  SEE  REMARKS.  DGEA2260 

IER  =  134,  IMPLIES  THAT  AFTER  SOME  INITIAL  DGEA2270 

SUCCESS,  THE  INTEGRATION  WAS  HALTED  EITHERDGEA22  80 

BY  REPEATED  ERROR  TEST  FAILURES  OR  BY  DGEA229  0 

A  TEST  ON  TOL.  SEE  REMARKS.  DGEA23  00 

IER  =  135,  IMPLIES  THAT  ONE  OF  THE  INPUT  DGEA2310 

PARAMETERS  N,X, H,XEND, TOL, METH, MITER,  OR  DGEA232  0 

INDEX  WAS  SPECIFIED  INCORRECTLY.  DGEA2330 

IER  =  13  6,  IMPLIES  THAT  INDEX  HAD  A  VALUE  DGEA234  0 

OF  -1  ON  INPUT,  BUT  THE  DESIRED  CHANGES  DGEA2350 

OF  PARAMETERS  WERE  NOT  IMPLEMENTED  DGEA23  60 

BECAUSE  XEND  WAS  NOT  BEYOND  X.  DGEA23  70 

INTERPOLATION  TO  X  =  XEND  WAS  PERFORMED.  DGEA2380 

TO  TRY  AGAIN,  SIMPLY  CALL  AGAIN  WITH  DGEA239  0 

INDEX  EQUAL  TO  -1  AND  A  NEW  VALUE  FOR  DGEA24  00 

XEND.  DGEA2410 

#  of  integration  steps  taken  + 

DGEA2420 

SINGLE  AND  DOUBLE /H3 2  DGEA2430 

SINGLE/H36,H48,H60  DGEA2440 


DGEA24  5  0 
REQD.  IMS L  ROUTINES  -  DGRCS , DGRIN, DGRPS , DGRST, LUDATF, LUELMF, LEQT1B,  DGEA2460 
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C  UERTST, UGETIO  DGEA24  70 

C  DGEA24  80 

C    NOTATION  -  INFORMATION  ON  SPECIAL  NOTATION  AND  DGEA24  9  0 

C  CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL  DGEA2500 

C  INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP  DGEA2510 

C  DGEA2520 

C   REMARKS   1.   THE  EXTERNAL  SUBROUTINE  FCNJ  IS  USED  ONLY  WHEN  DGEA2530 

C  INPUT  PARAMETER  MITER  IS  EQUAL  TO  1  OR  -1.  OTHERWISE,  DGEA2540 

C  A  DUMMY  FUNCTION  CAN  BE  USED.  THE  DUMMY  SUBROUTINE  DGEA255  0 

C  SHOULD  BE  OF  THE  FOLLOWING  FORM  DGEA2560 

C  SUBROUTINE  FCNJ  (N,X,Y,PD)  DGEA2570 

C  INTEGER  N  DGEA2580 

C  REAL  Y(N) ,PD(N,N) ,X  DGEA2590 

C  RETURN  DGEA2  600 

C  END  DGEA2610 

C  2.   AFTER  THE  INITIAL  CALL,  IF  A  NORMAL  RETURN  OCCURRED  DGEA2620 

C  (IER=0)  AND  A  NORMAL  CONTINUATION  IS  DESIRED,  SIMPLY  DGEA2  63  0 

C  RESET  XEND  AND  CALL  DGEAR  AGAIN.  ALL  OTHER  DGEA264  0 

C  PARAMETERS  WILL  BE  READY  FOR  THE  NEXT  CALL.  A  CHANGE  DGEA2650 

C  OF  PARAMETERS  WITH  INDEX  EQUAL  TO  -1  CAN  BE  MADE  DGEA2660 

C  AFTER  EITHER  A  SUCCESSFUL  OR  AN  UNSUCCESSFUL  RETURN.  DGEA2  670 

C  3.   THE  COMMON  BLOCKS  /DBAND/  AND  /GEAR/  NEED  TO  BE  DGEA2680 

C  PRESERVED  BETWEEN  CALLS  TO  DGEAR.  IF  IT  IS  NECESSARY  DGEA2690 

C  FOR  THE  COMMON  BLOCKS  TO  EXIST  IN  THE  CALLING  PROGRAM  DGEA2700 

C  THE  FOLLOWING  STATEMENTS  SHOULD  BE  INCLUDED  DGEA2710 

C  COMMON   /DBAND/  NLC,NUC  DGEA2  720 

C  COMMON   /GEAR/  DUMMY(48) ,SDUMMY(4) ,IDUMMY(38)  DGEA2730 

C  WHERE  DUMMY,  SDUMMY,  AND  IDUMMY  ARE  VARIABLE  NAMES  NOT  DGEA274  0 

C  USED  ELSEWHERE  IN  THE  CALLING  PROGRAM.   (FOR  DOUBLE  DGEA2750 

C  PRECISION  DUMMY  IS  TYPE  DOUBLE  AND  SDUMMY  IS  TYPE  REAL ) DGEA2 7 6 0 

C  4.   THE  CHOICE  OF  VALUES  FOR  METH  AND  MITER  MAY  REQUIRE  DGEA2770 

C  SOME  EXPERIMENTATION,  AND  ALSO  SOME  CONSIDERATION  OF  DGEA2  780 

C  THE  NATURE  OF  THE  PROBLEM  AND  OF  STORAGE  REQUIREMENTS.  DGEA2790 

C  THE  PRIME  CONSIDERATION  IS  STIFFNESS.  IF  DGEA2800 

C  THE  PROBLEM  IS  NOT  STIFF,  THE  BEST  CHOICE  IS  PROBABLY  DGEA2  810 

C  METH  =  1  WITH  MITER  =  0 .  IF  THE  PROBLEM  IS  STIFF  TO  A  DGEA2820 

C  SIGNIFICANT  DEGREE,  THEN  METH  SHOULD  BE  2  AND  MITER  DGEA2830 

C  SHOULD  BE  1,2,-1,-2  OR  3.  IF  THE  USER  HAS  NO  KNOWLEDGE  DGEA2840 

C  OF  THE  INHERENT  TIME  CONSTANTS  OF  THE  PROBLEM,  WITH  DGEA2850 

C  WHICH  TO  PREDICT  ITS  STIFFNESS,  ONE  WAY  TO  DETERMINE  DGEA2860 

C  THIS  IS  TO  TRY  METH  =  1  AND  MITER  =  0  FIRST,  AND  LOOK  DGEA2870 

C  AT  THE  BEHAVIOR  OF  THE  SOLUTION  COMPUTED  AND  THE  STEP  DGEA2  880 

C  SIZES  USED.  IF  THE  TYPICAL  VALUES  OF  H  ARE  MUCH  DGEA2890 

C  SMALLER  THAN  THE  SOLUTION  BEHAVIOR  WOULD  SEEM  TO  DGEA29  00 

C  REQUIRE  (THAT  IS,  MORE  THAN  100  STEPS  ARE  TAKEN  OVER  DGEA2910 

C  AN  INTERVAL  IN  WHICH  THE  SOLUTIONS  CHANGE  BY  LESS  DGEA2920 

C  THAN  ONE  PERCENT) ,  THEN  THE  PROBLEM  IS  PROBABLY  STIFF  DGEA2930 

C  AND  THE  DEGREE  OF  STIFFNESS  CAN  BE  ESTIMATED  FROM  THE  DGEA294  0 

C  VALUES  OF  H  USED  AND  THE  SMOOTHNESS  OF  THE  SOLUTION.  DGEA2950 

C  IF  THE  DEGREE  OF  STIFFNESS  IS  ONLY  SLIGHT,  IT  MAY  BE  DGEA2960 

C  THAT  METH=1  IS  MORE  EFFICIENT  THAN  METH=2 .  DGEA2970 

C  EXPERIMENTATION  WOULD  BE  REQUIRED  TO  DETERMINE  THIS.  DGEA2980 

C  REGARDLESS  OF  METH,  THE  LEAST  EFFECTIVE  VALUE  OF  DGEA2990 

C  MITER  IS  0,  AND  THE  MOST  EFFECTIVE  IS  1 , - 1 , 2 , OR  -2.  DGEA3000 

C  MITER  =  3  IS  GENERALLY  SOMEWHERE  IN  BETWEEN.  SINCE  DGEA3010 

C  THE  STORAGE  REQUIREMENTS  GO  UP  IN  THE  SAME  ORDER  AS  DGEA3  02  0 

C  EFFECTIVENESS,  TRADE-OFF  CONSIDERATIONS  ARE  DGEA3  03  0 

C  NECESSARY.  FOR  REASONS  OF  ACCURACY  AND  SPEED,  THE  DGEA3  04  0 

C  CHOICE  OF  ABS (MITER) =1  IS  GENERALLY  PREFERRED  TO  DGEA3050 

C  ABS (MITER) =2 ,  UNLESS  THE  SYSTEM  IS  FAIRLY  COMPLICATED  DGEA3060 

C  (AND  FCNJ  IS  THUS  NOT  FEASIBLE  TO  CODE)  .  THE  DGEA3  070 

C  ACCURACY  OF  THE  FCNJ  CALCULATION  CAN  BE  CHECKED  BY  DGEA3  080 
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COMPARISON  OF  THE  JACOB  IAN  WITH  THAT  GENERATED  WITH 
ABS (MITER) =2 .  IF  THE  JACOBIAN  MATRIX  IS  SIGNIFICANTLY 
DIAGONALLY  DOMINANT,  THEN  THE  OPTION  MITER  =  3  IS 
LIKELY  TO  BE  NEARLY  AS  EFFECTIVE  AS  ABS  (MITER)  =1  OR  2, 
AND  WILL  SAVE  CONSIDERABLE  STORAGE  AND  RUN  TIME. 
IT  IS  POSSIBLE,  AND  POTENTIALLY  QUITE  DESIRABLE,  TO 
USE  DIFFERENT  VALUES  OF  METH  AND  MITER  IN  DIFFERENT 
SUBINTERVALS  OF  THE  PROBLEM.  FOR  EXAMPLE,  IF  THE 
PROBLEM  IS  NON- STIFF  INITIALLY  AND  STIFF  LATER, 
METH  =  1  AND  MITER  =  0  MIGHT  BE  SET  INITIALLY,  AND 

METH  =  2  AND  MITER  =  1  LATER. 

THE  INITIAL  VALUE  OF  THE  STEP  SIZE,  H,  SHOULD  BE 
CHOSEN  CONSIDERABLY  SMALLER  THAN  THE  AVERAGE  VALUE 
EXPECTED  FOR  THE  PROBLEM,  AS  THE  FIRST-ORDER  METHOD 
WITH  WHICH  DGEAR  BEGINS  IS  NOT  GENERALLY  THE  MOST 
EFFICIENT  ONE.  HOWEVER,  FOR  THE  FIRST  STEP,  AS  FOR 
EVERY  STEP,  DGEAR  TESTS  FOR  THE  POSSIBILITY  THAT 
THE  STEP  SIZE  WAS  TOO  LARGE  TO  PASS  THE  ERROR  TEST 

(BASED  ON  TOL) ,  AND  IF  SO  ADJUSTS  THE  STEP  SIZE 
DOWN  AUTOMATICALLY.  THIS  DOWNWARD  ADJUSTMENT,  IF 
ANY,  IS  NOTED  BY  IER  HAVING  THE  VALUES  66  OR  67, 
AND  SUBSEQUENT  RUNS  ON  THE  SAME  OR  SIMILAR  PROBLEM 
SHOULD  BE  STARTED  WITH  AN  APPROPRIATELY  SMALLER 

VALUE  OF  H. 

SOME  OF  THE  VALUES  OF  INTEREST  LOCATED  IN  THE 

COMMON  BLOCK  /GEAR/  ARE 

A.  HUSED,  THE  STEP  SIZE  H  LAST  USED  SUCCESSFULLY 

(DUMMY  (8)  ) 

B.  NQUSED,  THE  ORDER  LAST  USED  SUCCESSFULLY 

( IDUMMY ( 6 ) ) 

C.  NSTEP,  THE  CUMULATIVE  NUMBER  OF  STEPS  TAKEN 

(IDUMMY(7) ) 

D.  NFE,  THE  CUMULATIVE  NUMBER  OF  FCN  EVALUATIONS 

(IDUMMY(8)  ) 

E.  NJE,  THE  CUMULATIVE  NUMBER  OF  JACOBIAN 
EVALUATIONS,  AND  HENCE  ALSO  OF  MATRIX  LU 
DECOMPOSITIONS  (IDUMMY (9)) 

THE  NORMAL  USAGE  OF  DGEAR  MAY  BE  SUMMARIZED  AS  FOLLOWS 

A.  SET  THE  INITIAL  VALUES  IN  Y. 

B.  SET  N,  X,  H,  TOL,  METH,  AND  MITER. 

C.  SET  XEND  TO  THE  FIRST  OUTPUT  POINT,  AND  INDEX  TO  1. 

D.  CALL  DGEAR 

E.  EXIT  IF  IER  IS  GREATER  THAN  128. 

F.  OTHERWISE,  DO  DESIRED  OUTPUT  OF  Y. 

G.  EXIT  IF  THE  PROBLEM  IS  FINISHED. 

H.  OTHERWISE,  RESET  XEND  TO  THE  NEXT  OUTPUT  POINT,  AND 

RETURN  TO  STEP  D. 
THE  ERROR  WHICH  IS  CONTROLLED  BY  WAY  OF  THE  PARAMETER 
TOL  IS  AN  ESTIMATE  OF  THE  LOCAL  TRUNCATION  ERROR,  THAT 
IS,  THE  ERROR  COMMITTED  ON  TAKING  A  SINGLE  STEP  WITH 
THE  METHOD,  STARTING  WITH  DATA  REGARDED  AS  EXACT.  THIS 
IS  TO  BE  DISTINGUISHED  FROM  THE  GLOBAL  TRUNCATION 
ERROR,  WHICH  IS  THE  ERROR  IN  ANY  GIVEN  COMPUTED  VALUE 
OF  Y(X)  AS  A  RESULT  OF  THE  LOCAL  TRUNCATION  ERRORS 
FROM  ALL  STEPS  TAKEN  TO  OBTAIN  Y(X) .  THE  LATTER  ERROR 
ACCUMULATES  IN  A  NON -TRIVIAL  WAY  FROM  THE  LOCAL 
ERRORS,  AND  IS  NEITHER  ESTIMATED  NOR  CONTROLLED  BY 
THE  ROUTINE.  SINCE  IT  IS  USUALLY  THE  GLOBAL  ERROR  THAT 
A  USER  WANTS  TO  HAVE  UNDER  CONTROL,  SOME 
EXPERIMENTATION  MAY  BE  NECESSARY  TO  GET  THE  RIGHT 
VALUE  OF  TOL  TO  ACHIEVE  THE  USERS  NEEDS.  IF  THE 
PROBLEM  IS  MATHEMATICALLY  STABLE,  AND  THE  METHOD  USED 


DGEA3090 

DGEA3100 

DGEA3110 

DGEA3120 

DGEA313  0 

DGEA3140 

DGEA3150 

DGEA3160 

DGEA3170 

DGEA3180 

DGEA3190 

DGEA3200 

DGEA3210 

DGEA3220 

DGEA3230 

DGEA3240 

DGEA3250 

DGEA3260 

DGEA3270 

DGEA3280 

DGEA3290 

DGEA3300 

DGEA3310 

DGEA3320 

DGEA3330 

DGEA3340 

DGEA3350 

DGEA3360 

DGEA3370 

DGEA3380 

DGEA3390 

DGEA34  00 

DGEA3410 

DGEA3420 

DGEA3430 

DGEA344  0 

DGEA3450 

DGEA34  60 

DGEA34  70 

DGEA34  80 

DGEA3490 

DGEA3500 

DGEA3510 

DGEA352  0 

DGEA3530 

DGEA354  0 

DGEA3550 

DGEA3560 

DGEA3570 

DGEA3580 

DGEA359  0 

DGEA3  600 

DGEA3  610 

DGEA3  620 

DGEA3630 

DGEA3  64  0 

DGEA3  650 

DGEA3660 

DGEA3670 

DGEA3680 

DGEA369  0 

DGEA3700 
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IS  APPROPRIATELY  STABLE,  THEN  THE  GLOBAL  ERROR  AT  A  DGEA3  710 

GIVEN  X  SHOULD  VARY  SMOOTHLY  WITH  TOL  IN  A  MONOTONE  DGEA3720 

INCREASING  MANNER.  DGEA3  73  0 

IF  THE  ROUTINE  RETURNS  WITH  IER  VALUES  OF  132,  133,  ~GEA3740 

OR  134,  THE  USER  SHOULD  CHECK  TO  SEE  IF  TOO  MUCH  DGEA3750 

ACCURACY  IS  BEING  REQUIRED.  THE  USER  MAY  WISH  TO  DGEA3760 

SET  TOL  TO  A  LARGER  VALUE  AND  CONTINUE.  ANOTHER  DGEA377  0 

POSSIBLE  CAUSE  OF  THESE  ERROR  CONDITIONS  IS  AN  DGEA3  780 

ERROR  IN  THE  CODING  OF  THE  EXTERNAL  FUNCTIONS  FCN  DGEA379  0 

OR  FCNJ.  IF  NO  ERRORS  ARE  FOUND,  IT  MAY  BE  NECESSARY  DGEA3800 

TO  MONITOR  INTERMEDIATE  QUANTITIES  GENERATED  BY  THE  DGEA3  810 
ROUTINE.  THESE  QUANTITIES  ARE  STORED  IN  THE  WORK  VECTORDGEA3  82  0 
WK  AND  INDEXED  BY  SPECIFIC  ELEMENTS  IN  THE  COMMON  BLOCKDGEA3  83  0 

/GEAR/.  IF  IER  IS  132  OR  134,  THE  COMPONENTS  CAUSING  DGEA3840 

THE  ERROR  TEST  FAILURE  CAN  BE  IDENTIFIED  FROM  LARGE  DGEA3850 

VALUES  OF  THE  QUANTITY  DGEA3  860 

WK(IDUMMY(11)+I) /WK(I) ,  FOR  1=1 , . . . , N .  DGEA3870 

ONE  CAUSE  OF  THIS  MAY  BE  A  VERY  SMALL  BUT  NONZERO  DGEA3  880 

INITIAL  VALUE  OF  ABS(Y(I)) .  DGEA3  89  0 

IF  IER  IS  133,  SEVERAL  POSSIBILITIES  EXIST.  DGEA3900 
IT  MAY  BE  INSTRUCTIVE  TO  TRY  DIFFERENT  VALUES  OF  MITER .DGEA39 10 


COPYRIGHT 


WARRANTY 


ALTERNATIVELY,  THE  USER  MIGHT  MONITOR  SUCCESSIVE 
CORRECTOR  ITERATES  CONTAINED  IN  WK ( IDUMMY ( 12 ) +1 ) ,  FOR 
1=1, . . . ,N.  ANOTHER  POSSIBILITY  MIGHT  BE  TO  MONITOR 
THE  JACOBIAN  MATRIX,  IF  ONE  IS  USED,  STORED,  BY 
COLUMN,  IN  WK(IDUMMY(10)+I) ,  FOR  1=1 , . . . , N*N  IF 
ABS (MITER)  IS  EQUAL  TO  1  OR  2 ,  OR  FOR  1=1, . . .,N  IF 
MITER  IS  EQUAL  TO  3 . 

-  1984  BY  IMSL,  INC.   ALL  RIGHTS  RESERVED. 


EXTERNAL 
COMMON  /DBAND/ 
COMMON  /GEAR/ 


DGEA3920 

DGEA393  0 

DGEA394  0 

DGEA3950 

DGEA39  60 

DGEA3970 

DGEA39  80 

DGEA3990 

DGEA4000 

DGEA4  010 

IMSL  WARRANTS  ONLY  THAT  IMSL  TESTING  HAS  BEEN  DGEA4  020 

APPLIED  TO  THIS  CODE.   NO  OTHER  WARRANTY,    DGEA4030 

EXPRESSED  OR  IMPLIED,  IS  APPLICABLE.         DGEA4  04  0 

DGEA4  050 

DGEA4060 

DGEA4  07  0 
DGEA4  080 

+ 
DGEA4100 

+ 
DGEA4120 
DGEA413  0 
DGEA414  0 
,DGEA4150 
DGEA4160 
DGEA417  0 
DGEA4180 
DGEA419  0 
DGEA42  00 
DGEA4210 
DGEA422  0 
DGEA423  0 
DGEA4  24  0 
DGEA4  250 
DGEA4260 
DGEA4270 
DGEA42  80 
DGEA429  0 
DGEA4300 
DGEA4  310 
DGEA432  0 


SUBROUTINE  DGEAR   ( N , FCN , FCNJ , X , H , Y , XEND , TOL , METH , MITER , INDEX , 
1  IWK,WK, IER, step) 

SPECIFICATIONS  FOR  ARGUMENTS 
INTEGER  N,METH, MITER, INDEX, IWK(l) , IER, step 

DOUBLE  PRECISION   X, H, Y (N) , XEND, TOL, WK (1) 

SPECIFICATIONS  FOR  LOCAL  VARIABLES 
INTEGER  NERROR, NSAVE1 , NSAVE2 , NPW, NY, NC , MFC , KFLAG , 

1  JSTART , NSQ , NQUSED , NSTEP , NFE , NJE , I , NO , NHCUT , KGO 

2  JER,KER,NN,NEQUIL, IDUMMY(21) , NLC , NUC 
REAL  SDUMMY(4) 
DOUBLE  PRECISION   T, HH, HMIN, HMAX, EPSC,UROUND, EPSJ, HUSED, TOUTP, 

1  AYI,D,DN,SEPS, DUMMY (39) 

FCN, FCNJ 
NLC, NUC 
T , HH , HMIN , HMAX , EPSC , UROUND , EPSJ , HUSED , DUMMY , 

1  TOUTP , SDUMMY , NC , MFC , KFLAG , JSTART , NSQ , NQUSED , 

2  NSTEP , NFE , NJE , NPW, NERROR, NSAVE1 , NSAVE2 , NEQUIL, 

3  NY, IDUMMY, NO, NHCUT 
DATA  SEPS/Z3410000000000000/ 

FIRST  EXECUTABLE  STATEMENT 
IF  (MITER. GE.O)  NLC  =  -1 
KER  =  0 
JER  =  0 
URCUND  =  SEPS 

COMPUTE  WORK  VECTOR  INDICIES 
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NERROR  =  N 
NSAVE1  =  NERROR+N 
NSAVE2  =  NSAVE1+N 
NY  =  NSAVE2+N 
IF  (METH.EQ.l) 
IF  (METH.EQ.2) 
NPW  =  NEQUIL  + 
IF  (MITER. EQ.O 


NEQUIL  =  N5T+13*N 

NEQUIL  =  NY+6*N 

N 

OR. MITER. EQ. 3)  NPW  =  NEQUIL 


MFC  =  10*METH+IABS (MITER) 


10 


15 


20 


CHECK  FOR  INCORRECT  INPUT  PARAMETERS 


IF 
IF 
IF 
IF 
IF 
IF 
IF 
IF 
IF 
IF 


(MITER. LT. -2. OR. MITER. GT. 3)  GOTO  85 

(METH.NE.1.AND.METH.NE.2)  GO  TO  85 

(TOL.LE.0.D0)  GO  TO  85 

(N.LE.O)  GO  TO  85 

( (X-XEND)*H.GE.0.D0)  GO  TO  85 

(INDEX. EQ.O)  GO  TO  10 

(INDEX.EQ.2)  GO  TO  15 


( INDEX. EQ 
( INDEX. EQ 
( INDEX. NE 


-1)  GO  TO  20 
3)  GO  TO  25 
1)  GO  TO  85 


IF  INITIAL  VALUES  OF  YMAX  OTHER  THAN 
THOSE  SET  BELOW  ARE  DESIRED,  THEY 
SHOULD  BE  SET  HERE.  ALL  YMAX (I) 
MUST  BE  POSITIVE.  IF  VALUES  FOR 
HMIN  OR  HMAX,  THE  BOUNDS  ON 
DABS (HH) ,  OTHER  THAN  THOSE  BELOW 
ARE  DESIRED,  THEY  SHOULD  BE  SET 
BELOW . 


DO 


5  1=1, N 

WK(I)  =  DABS(Y(I)  ) 

IF  (WK(I) .EQ.O. DO)  WK(I)  =  1.D0 

WK(NY+I)  =  Y(I) 
CONTINUE 
NC  =  N 
T  =  X 
HH  =  H 

IF  ( (T+HH) .EQ.T)  KER  = 
HMIN  =  DABS (H) 
HMAX  =  DABS (X-XEND) *10 
EPSC  =  TOL 
JSTART  =  0 
NO  =  N 
NSQ  =  N0*N0 
EPSJ  =  DSQRT(UROUND) 
NHCUT  =  0 
DUMMY (2)  =  1.0D0 
DUMMY (14)  =  1.0D0 
GO  TO  30 


33 


.DO 


HMAX  =  DABS(XEND-TOUTP) *10 
GO  TO  45 


TOUTP  IS  THE  PREVIOUS  VALUE  OF  XEND 

FOR  USE  IN  HMAX. 
DO 


HMAX  =  DABS (XEND -TOUTP) *10. DO 
IF  ( (T-XEND) *HH.GE.0.D0)  GOTO 
GO  TO  50 


95 


IF  ( (T-XEND) *HH.GE 
JSTART  =  -1 
NC  =  N 
EPSC  =  TOL 


0.D0)  GO  TO  90 


DGEA4330 
DGEA4340 
DGEA4350 
DGEA4360 
DGEA4  370 
DGEA4380 
DGEA4390 
DGEA44  00 
DGEA4410 
DGEA4420 
DGEA4430 
DGEA4440 
DGEA4450 
DGEA44  60 
DGEA4470 
DGEA44  80 
DGEA4490 
DGEA4500 
DGEA4510 
DGEA4520 
DGEA4  530 
DGEA454  0 
DGEA4  550 
DGEA4560 
DGEA4570 
DGEA4  580 
DGEA4590 
DGEA4  600 
DGEA4610 
DGEA4620 
DGEA4  630 
DGEA4  64  0 
DGEA4  650 
DGEA4  660 
DGEA4  670 
DGEA4  680 
DGEA4  690 
DGEA4  700 
DGEA4  710 
DGEA4  720 
DGEA4  73  0 
DGEA4  74  0 
DGEA4  750 
DGEA4  760 
DGEA4  770 
DGEA4  780 
DGEA4  79  0 
DGEA4800 
DGEA4  810 
DGEA4  82  0 
DGEA4  830 
DGEA4  84  0 
DGEA4  85  0 
DGEA4  860 
DGEA4  870 
DGEA4  880 
DGEA4  890 
DGEA49  00 
DGEA4910 
DGEA4920 
DGEA4930 
DGEA494  0 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


25  IF  ( (T+HH) .EQ.T)  KER  =  33 

write  (*,*),' error  code  =  '  ,  Jeer 

30  NN  =  NO 

step  =  step  +  1 

write (*,*)' step  =  ' , step 

CALL  DGRST  (FCN, FCNJ, WK (NY+1) , WK, WK (NERROR+1) , WK (NSAVE1+1) , 
1  WK(NSAVE2+1) ,WK(NPW+1) , WK (NEQUIL+1) , IWK, NN, step) 


KGO  =  1-KFLAG 

GO  TO  (35,55,70,80) , 

3  5  CONTINUE 


40 


45 


50 


KGO 


KFLAG  =0,  -1,  -2,  -3 

NORMAL  RETURN  FROM  INTEGRATOR.  THE 
WEIGHTS  YMAX(I)  ARE  UPDATED.  IF 
DIFFERENT  VALUES  ARE  DESIRED,  THEY 
SHOULD  BE  SET  HERE.  A  TEST  IS  MADE 
FOR  TOL  BEING  TOO  SMALL  FOR  THE 
MACHINE  PRECISION.  ANY  OTHER  TESTS 
OR  CALCULATIONS  THAT  ARE  REQUIRED 
AFTER  EVERY  STEP  SHOULD  BE 
INSERTED  HERE.  IF  INDEX  =3,  Y  IS 
SET  TO  THE  CURRENT  SOLUTION  ON 
RETURN.  IF  INDEX  =2,  HH  IS 
CONTROLLED  TO  HIT  XEND  (WITHIN 
ROUNDOFF  ERROR) ,  AND  THEN  THE 
CURRENT  SOLUTION  IS  PUT  IN  Y  ON 
RETURN.  FOR  ANY  OTHER  VALUE  OF 
INDEX,  CONTROL  RETURNS  TO  THE 
INTEGRATOR  UNLESS  XEND  HAS  BEEN 
REACHED.  THEN  INTERPOLATED  VALUES 
OF  THE  SOLUTION  ARE  COMPUTED  AND 
STORED  IN  Y  ON  RETURN. 
IF  INTERPOLATION  IS  NOT 
DESIRED,  THE  CALL  TO  DGRIN  SHOULD 
BE  REMOVED  AND  CONTROL  TRANSFERRED 
TO  STATEMENT  95  INSTEAD  OF  105. 

D  =  0.D0 

DO  40  1  =  1, N 

AYI  =  DABS (WK(NY+I) ) 
WK(I)  =  DMAX1 (WK(I) ,AYI) 

D  =  D+ (AYI/WK(I) ) **2 

D  =  D* (UROUND/TOL) **2 

DN  =  N 

IF  (D.GT.DN)  GO  TO  75 

IF  (INDEX.EQ.3)  GO  TO  95 

IF  (INDEX.EQ.2)  GO  TO  50 

IF  ( (T-XEND) *HH.LT.0.D0)  GO  TO  25 

NN  =  NO 

CALL  DGRIN  (XEND , WK (NY+1) , NN, Y) 

X  =  XEND 

GO  TO  105 

IF  (( (T+HH) -XEND) *HH.LE.0. DO)  GO  TO  25 

IF  (DABS (T-XEND) .LE.UROUND*DMAXl (10. D0*DABS (T) ,HMAX) )  GO  TO  95 

IF  ( (T-XEND) *HH.GE.0. DO)  GO  TO  95 

HH  =  (XEND-T) * (1.D0-4 .D0*UROUND) 

JSTART  =  -1 

GO  TO  25 

ON  AN  ERROR  RETURN  FROM  INTEGRATOR, 
AN  IMMEDIATE  RETURN  OCCURS  IF 
KFLAG  =  -2,  AND  RECOVERY  ATTEMPTS 


DGEA4950 
DGEA4  9  60 

+ 
DGEA4  970 
DGEA49  8  0 

+ 

+ 
DGEA4  99  0 

+ 
DGEA5010 
DGEA5020 
DGEA503  0 
DGEA5  04  0 
DGEA5050 
DGEA5060 
DGEA5070 
DGEA5080 
DGEA5090 
DGEA5100 
DGEA5110 
DGEA5120 
DGEA513  0 
DGEA514  0 
DGEA5150 
DGEA5160 
DGEA5170 
DGEA5180 
DGEA5190 
DGEA5200 
DGEA5210 
DGEA5220 
DGEA5230 
DGEA5240 
DGEA5250 
DGEA5260 
DGEA5270 
DGEA52  80 
DGEA5290 
DGEA5300 
DGEA5310 
DGEA5320 
DGEA53  3  0 
DGEA534  0 
DGEA53  50 
DGEA53  60 
DGEA5370 
DGEA53  80 
DGEA539  0 
DGEA54  00 
DGEA5410 
DGEA542  0 
DGEA54  3  0 
DGEA544  0 
DGEA54  50 
DGEA54  60 
DGEA54  70 
DGEA54  80 
DGEA54  9  0 
DGEA5500 
DGEA5510 
DGEA552  0 
DGEA553  0 
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c 
c 
c 
c 


55 
60 


65 


JER  =  66 

IF  (NHCUT.EQ.10)  GO  TO  65 

NHCUT  =  NHCUT+1 

HMIN  =  HMIN*.1D0 

HH  =  HH* . 1D0 

JSTART  =  -1 

GO  TO  25 

IF  (JER.EQ.66)  JER  =  132 
IF  (JER.EQ.67)  JER  =  133 
GO  TO  95 


ARE  MADE  OTHERWISE.  TO  RECOVER   HH 
AND  HMIN  ARE  REDUCED  BY  A  FACTOR 
OF  .1  UP  TO  10  TIMES  BEFORE  GIVING 


70  JER  =  134 
GO  TO  95 

JER  =  134 
KFLAG  =  -2 
GO  TO  95 

JER  =67 
GO  TO  60 

JER  =135 
GO  TO  110 

JER  =  136 
NN  =  NO 

CALL  DGRIN  (XEND, WK (NY+1)  NN  Y) 
X  =  XEND  '   ' 

GO  TO  110 


95 


75 


80 


85 


90 


INDEX  =  KFLAG 


X  =  T 

DO  100  1=1, N 
Y(I)  =  WK(NY+I) 
IF  (JER. LT. 128) 
TOUTP  =  X 

IF  (KFLAG. EQ.0)  H  =  HUSED 
IF  (KFLAG. NE.0)  H  =  HH 
IER  =  MAX0 (KER, JER) 
CONTINUE 

IF  (KER. NE.0. AND.JER.LT.  128)  CALL  UERTST  (KER  fiHnPPAP  i 
9005  REFTUR^R-NE-0)  «"■*»«*  (OER.6HDGeSST  ^'^^    > 
END 


100 
105 


110 
9000 


DGEA554  0 

DGEA5550 

DGEA5560 

DGEA5570 

DGEA5580 

DGEA5590 

DGEA5600 

DGEA5610 

DGEA5  62  0 

DGEA5  630 

DGEA564  0 

DGEA5650 

DGEA5660 

DGEA5670 

DGEA5680 

DGEA5690 

DGEA5700 

DGEA5710 

DGEA5720 

DGEA573  0 

DGEA574  0 

DGEA5750 

DGEA5760 

DGEA5770 

DGEA5780 

DGEA5790 

DGEA5800 

DGEA5810 

DGEA5820 

DGEA5830 

DGEA5840 

DGEA5850 

DGEA5860 

DGEA5870 

DGEA5880 

DGEA5890 

DGEA5900 

DGEA5910 

DGEA5920 

DGEA5930 

DGEA5940 

DGEA5950 

DGEA59  60 

DGEA5970 

DGEA5980 

DGEA5990 

DGEA6  000 

DGEA6010 
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IMSL  ROUTINE  NAME    -  DGRST  DGRS0010 

DGRS0020 
modified  to  print  sheath  and  diagnostic  output  to  files  "sheatha.dat"  + 
and  "diag.dat" 

COMPUTER 


IBM/DOUBLE 


LATEST  REVISION 


PURPOSE 


PRECISION/HARDWARE 


REQD.  IMSL  ROUTINES 


NOTATION 


COPYRIGHT 


WARRANTY 


•  JUNE  1,  1982 

NUCLEUS  CALLED  ONLY  BY  IMSL  SUBROUTINE  DGEAR 

-  SINGLE  AND  DOUBLE /H3 2 
SINGLE /H3 6 , H4  8 , H6  0 

DGRCS , DGRPS , LUDATF , LUELMF , LEQT1B , UERTST , 
UGETIO 

INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 

19  82  BY  IMSL,  INC.  ALL  RIGHTS  RESERVED. 


SUBROUTINE  DGRST 


INTEGER 

DOUBLE  PRECISION 


INTEGER 


REAL 

DOUBLE  PRECISION 


EXTERNAL 
COMMON  /DBAND/ 
COMMON  /GEAR/ 


open (unit=8 , file 
open (unit=9 , file 
KFLAG  =  0 
TO'iD  =  T 


DGRS0050 
DGRS0060 
DGRS0070 
DGRS0080 
DGRS0090 
DGRS0100 
DGRS0110 
DGRS0120 
DGRS0130 
DGRS0140 
DGRS0150 
DGRS0160 
DGRS0170 
DGRS0180 
DGRS0190 
DGRS0200 
DGRS0210 
DGRS0220 
IMSL  WARRANTS  ONLY  THAT  IMSL  TESTING  HAS  BEEN  DGRS023  0 
APPLIED  TO  THIS  CODE.  NO  OTHER  WARRANTY,  DGRS0240 
EXPRESSED  OR  IMPLIED,  IS  APPLICABLE.         DGRS0250 

DGRS0260 

DGRS0270 

DGRS0280 
DGRS0290 

+ 
DGRS0310 
DGRS0320 
DGRS0330 

+ 
DGRS0350 
DGRS0360 
DGRS0370 
DGRS0380 
DGRS0390 
DGRS0400 
DGRS0410 
DGRS0420 
DGRS0430 
,DGRS0440 

+ 
DGRS0460 
DGRS04  70 
DGRS0480 
DGRS0490 
DGRS0500 
DGRS0510 
DGRS0520 
DGRS0530 
DGRS0540 
DGRS0550 
DGRS0560 
DGRS0570 

+ 

+ 
DGRS0580 
DGRS0590 
DGRS0600 


( FCN , FCNJ , Y , YMAX , ERROR , SAVE 1 , SAVE2 , PW , EQUIL , 
IPIV,NO,step) 

SPECIFICATIONS  FOR  ARGUMENTS 
IPIV(l) ,N0 

Y (NO , 1 ) , YMAX ( 1 ) , ERROR ( 1 ) , SAVE1 ( 1 ) , SAVE2 ( 1 ) , 
PW ( 1 ) , EQUIL ( 1 ) , epr ime , epr ime ( 2 ) 

SPECIFICATIONS  FOR  LOCAL  VARIABLES 
N , MF , KFLAG , JSTART , NQUSED , NSTEP , NFE , NJE , NSQ , 
I , METH , MITER , NQ , L , IDOUB , MFOLD , NOLD , IRET ,  MEO , 
MIO , IWEVAL, MAXDER, LMAX , IREDO , J , NSTEP J , Jl , J2 , 
M, IER, NEWQ, NPW, NERROR, NSAVE1 , NSAVE2 , NEQUIL, NY, 
MITER1 , IDUMMY ( 2 ) , NLC , NUC , NWK , JER 
TQ(4) 

T,H,HMIN,HMAX,EPS,UROUND,HUSED,EL(13) , OLDLO , 
TOLD , RMAX , RC , CRATE , EPSOLD , HOLD , FN , EDN , E , EUP , 
BND , RH, Rl , CON, R, HLO , RO , D , PHLO , PR3 , Dl , ENQ3 , ENQ2 
PR2 , PR1 , ENQ1 , EPS J , DUMMY , t cum 
FCN, FCNJ 
NLC , NUC 

T , H , HMIN , HMAX , EPS , UROUND , EPS J , HUSED , 
EL , OLDLO , TOLD , RMAX , RC , CRATE , EPSOLD , HOLD , FN , 
EDN, E , EUP , BND , RH , Rl , R, HLO , RO , D , PHLO , PR3 , Dl , 
ENQ3 , ENQ2 , PR2 , PR1 , ENQ1 , DUMMY , TQ , 
N , MF , KFLAG , JSTART , NSQ , NQUSED , NSTEP , NFE , NJE , 
NPW, NERROR, NS AVE 1 , NSAVE2 , NEQUIL, NY, 
I , METH , MITER , NQ , L , IDOUB , MFOLD , NOLD , IRET , MEO , 
MIO , IWEVAL , MAXDER, LMAX , IREDO , J , NSTEP J , Jl , J2 , 
M, NEWQ, IDUMMY 

FIRST  EXECUTABLE  STATEMENT 
=' sheatha.dat' , status=' unknown' ) 
=' diag.dat' , status=' unknown' ) 


THIS  ROUTINE  PERFORMS  ONE  STEP  OF 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


IF  (JSTART.GT.O)  GO  TO  50 
IF  (JSTART.NE.O)  GO  TO  10 


THE  INTEGRATION  OF  AN  INITIAL 
VALUE  PROBLEM  FOR  A  SYSTEM  OF 
ORDINARY  DIFFERENTIAL  EQUATIONS. 


C 
C 
C 
C 
C 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


ON  THE  FIRST  CALL,  THE  ORDER  IS  SET 
TO  1  AND  THE  INITIAL  YDOT  IS 
CALCULATED.  RMAX  IS  THE  MAXIMUM 
RATIO  BY  WHICH  H  CAN  BE  INCREASED 
IN  A  SINGLE  STEP.  IT  IS  INITIALLY 
1.E4  TO  COMPENSATE  FOR  THE  SMALL 
INITIAL  H,  BUT  THEN  IS  NORMALLY 
EQUAL  TO  10.  IF  A  FAILURE  OCCURS 
(IN  CORRECTOR  CONVERGENCE  OR  ERROR 
TEST) ,  RMAX  IS  SET  AT  2  FOR  THE 
NEXT  INCREASE. 

CALL  FCN  (N,T,Y,SAVEl,eprime,eprime2) 

DO  5  1=1,  N 

Y(I, 2)  =  H*SAVE1 (I) 

METH  =  MF/10 

MITER  =  MF-10*METH 

NQ  =  1 

L  =  2 

IDOUB  =  3 

RMAX  =  1.D4 

RC  =  0.D0 

CRATE  =  1.D0 

HOLD  =  H 

MFOLD  =  MF 

NSTEP  =  0 

NSTEPJ  =  0 

NFE  =  1 

NJE  =  0 

IRET  =  3 

GO  TO  15 

IF  THE  CALLER  HAS  CHANGED  METH, 
DGRCS  IS  CALLED  TO  SET  THE 
COEFFICIENTS  OF  THE  METHOD.  IF  THE 
CALLER  HAS  CHANGED  N,  EPS,  OR 
METH,  THE  CONSTANTS  E 
AND  BND  MUST  BE  RESET 
COMPARISON  FOR  ERRORS 
CURRENT  ORDER  NQ .  EUP 
FOR  INCREASING  THE  ORDER, 
DECREASING  THE  ORDER.  BND 


EDN,  EUP, 
.  E  IS  A 
OF  THE 
IS  TO  TEST 
EDN  FOR 
IS  USED 

TO  TEST  FOR  CONVERGENCE  OF  THE 
CORRECTOR  ITERATES.  IF  THE  CALLER 
HAS  CHANGED  H,  Y  MUST  BE  RESCALED . 
IF  H  OR  METH  HAS  BEEN  CHANGED, 
IDOUB  IS  RESET  TO  L  +  1  TO  PREVENT 
FURTHER  CHANGES  IN  H  FOR  THAT  MANY 
STEPS. 


10  IF  (MF.EQ. MFOLD)  GO  TO  25 
MEO  =  METH 
MIO  =  MITER 
METH  =  MF/10 
MITER  =  MF-10*METH 
MFOLD  =  MF 

IF  (MITER. NE. MIO)  IWEVAL  =  MITER 
IF  (METH. EQ. MEO)  GO  TO  25 
IDOUB  =  L+l 
IRET  =  1 


DGRS0610 
DGRS0620 
DGRS0630 
DGRS0640 
DGRS0650 
DGRS0660 
DGRS0670 
DGRS0680 
DGRS0690 
DGRS0700 
DGRS0710 
DGRS0720 
DGRS0730 
DGRS0740 
DGRS0750 
DGRS0760 

+ 
DGRS0780 
DGRS0790 
DGRS0800 
DGRS0810 
DGRS0820 
DGRS083  0 
DGRS0840 
DGRS0850 
DGRS0860 
DGRS0870 
DGRS0880 
DGRS089  0 
DGRS0900 
DGRS0910 
DGRS0920 
DGRS0930 
DGRS0940 
DGRS0950 
DGRS09  60 
DGRS0970 
DGRS0980 
DGRS0990 
DGRS1000 
DGRS1010 
DGRS1020 
DGRS1030 
DGRS104  0 
DGRS1050 
DGRS1060 
DGRS1070 
DGRS1080 
DGRS1090 
DGRS1100 
DGRS1110 
DGRS1120 
DGRS1130 
DGRS114  0 
DGRS1150 
DGRS1160 
DGRS1170 
DGRS1180 
DGRS1190 
DGRS1200 
DGRS1210 
DGRS1220 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


c 
c 
c 
c 
c 
c 
c 
c 
c 


15  CALL  DGRCS  (METH, NQ, EL, TQ,MAXDER) 
LMAX  =  MAXDER+1 
RC  =  RC*EL(1) /OLDLO 
OLDLO  =  EL(1) 
20  FN  =  N 

EDN  =  FN* (TQ(1) *EPS) **2 
E  =  FN* (TQ(2) *EPS) **2 
EUP  =  FN*  (TQ(3)  *EPS)**2 
BND  =  FN* (TQ(4) *EPS) **2 
EPSOLD  =  EPS 
NOLD  =  N 

GO  TO  (3  0,35,50) ,     IRET 
25  IF  ( (EPS. EQ. EPSOLD) .AND. (N.EQ. NOLD) )  GOTO  30 
IF  (N.EQ. NOLD)  IWEVAL  =  MITER 
IRET  =  1 
GO  TO  2  0 
30  IF  (H.EQ.HOLD)  GO  TO  50 
RH  =  H/HOLD 
H  =  HOLD 
IREDO  =  3 
GO  TO  40 
35  RH  =  DMAX1 (RH, HMIN/DABS (H) ) 
40  RH  =  DMIN1 (RH,HMAX/DABS (H) , RMAX) 
Rl  =  1.D0 
DO  45  J=2,L 

Rl  =  R1*RH 
DO  45  1=1, N 
45  Y(I,  J)  =  Y(I,  J)  *R1 
H  =  H*RH 
RC  =  RC*RH 
IDOUB  =  L+l 
IF  ( IREDO. EQ.0)  GO  TO  285 

THIS  SECTION  COMPUTES  THE  PREDICTED 
VALUES  BY  EFFECTIVELY  MULTIPLYING 
THE  Y  ARRAY  BY  THE  PASCAL  TRIANGLE 
MATRIX.  RC  IS  THE  RATIO  OF  NEW  TO 
OLD  VALUES  OF  THE  COEFFICIENT 
H*EL(1) .  WHEN  RC  DIFFERS  FROM  1  BY 
MORE  THAN  3  0  PERCENT,  OR  THE 
CALLER  HAS  CHANGED  MITER,  IWEVAL 
IS  SET  TO  MITER  TO  FORCE  THE 
PARTIALS  TO  BE  UPDATED,  IF 
PARTIALS  ARE  USED.  IN  ANY  CASE, 
THE  PARTIALS  ARE  UPDATED  AT  LEAST 
EVERY  20 -TH  STEP. 
50  IF  (DABS(RC-l.DO) .GT.0.3D0)  IWEVAL  =  MITER 
IF  (NSTEP.GE.NSTEPJ+2  0)  IWEVAL  =  MITER 
T  =  T+H 
DO  55  J1=1,NQ 
DO  55  J2=J1,NQ 

J  =  (NQ+J1) -J2 
DO  55  1=1, N 
55  Y(I, J)  =  Y(I,  J)+Y(I,  J+l) 

UP  TO  3  CORRECTOR  ITERATIONS  ARE 
TAKEN.  A  CONVERGENCE  TEST  IS  MADE 
ON  THE  R.M.S.  NORM  OF  EACH 
CORRECTION,  USING  BND,  WHICH  IS 
DEPENDENT  ON  EPS.  THE  SUM  OF  THE 
CORRECTIONS  IS  ACCUMULATED  IN  THE 
VECTOR  ERROR (I) .  THE  Y  ARRAY  IS 
NOT  ALTERED  IN  THE  CORRECTOR  LOOP. 
THE  UPDATED  Y  VECTOR  IS  STORED 


DGRS1230 

DGRS1240 

DGRS1250 

DGRS1260 

DGRS1270 

DGRS1280 

DGRS1290 

DGRS1300 

DGRS1310 

DGRS1320 

DGRS133  0 

DGRS1340 

DGRS1350 

DGRS1360 

DGRS1370 

DGRS1380 

DGRS1390 

DGRS14  00 

DGRS1410 

DGRS1420 

DGRS14  3  0 

DGRS1440 

DGRS1450 

DGRS1460 

DGRS14  70 

DGRS14  80 

DGRS1490 

DGRS1500 

DGRS1510 

DGRS1520 

DGRS1530 

DGRS154  0 

DGRS1550 

DGRS1560 

DGRS1570 

DGRS1580 

DGRS1590 

DGRS1600 

DGRS1610 

DGRS1620 

DGRS1630 

DGRS1640 

DGRS1650 

DGRS1660 

DGRS1670 

DGRS1680 

DGRS1690 

DGRS1700 

DGRS1710 

DGRS1720 

DGRS1730 

DGRS174  0 

DGRS1750 

DGRS1760 

DGRS1770 

DGRS1780 

DGRS1790 

DGRS1800 

DGRS1810 

DGRS1820 

DGRS1830 

DGRS1840 
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C  TEMPORARILY  IN  SAVE1 .               DGRS1850 

60  DO  65  1=1, N  DGRS1860 

65  ERROR (I)  =  0.D0  DGRS1870 

M  =  0  DGRS1880 

CALL  FCN  (N,T, Y, SAVE2 , eprime , eprime2)  + 

NFE  =  NFE+1  DGRS1900 

IF  ( IWEVAL . LE . 0 )  GO  TO  95  DGRS1910 

C  IF  INDICATED,  THE  MATRIX  P  =  I  -      DGRS1920 

C  H*EL(1)*J  IS  REEVALUATED  BEFORE     DGRS1930 

C  STARTING  THE  CORRECTOR  ITERATION.   DGRS1940 

C  IWEVAL  IS  SET  TO  0  AS  AN  INDICATOR  DGRS1950 

C  THAT  THIS  HAS  BEEN  DONE.  IF  MITER   DGRS1960 

C  =  1  OR  2,  P  IS  COMPUTED  AND         DGRS1970 

C  PROCESSED  IN  PSET.  IF  MITER  =  3,    DGRS1980 

C  THE  MATRIX  USED  IS  P  =  I  -          DGRS199  0 

C  H*EL(1)*D,  WHERE  D  IS  A  DIAGONAL    DGRS2000 

C  MATRIX.                             DGRS2010 

IWEVAL  =  0  DGRS2020 

RC  =  1.D0  DGRS2  03  0 

NJE  =  NJE+1  DGRS204  0 

NSTEPJ  =  NSTEP  DGRS2050 

GOTO  (75,70,80),  MITER  DGRS2060 

7  0  NFE  =  NFE+N  DGRS2070 

75  CON  =  -H*EL(1)  DGRS2080 

MITER1  =  MITER  DGRS2  09  0 

CALL  DGRPS  (FCN, FCNJ, Y, NO , CON,MITERl , YMAX, SAVE1 , SAVE2 , PW, EQUIL,   DGRS2100 

1  IPIV,IER)  DGRS2110 

IF  (IER.NE.O)  GO  TO  155  DGRS2120 

GO  TO  125  DGRS2130 

80  R  =  EL(1)*.1D0  DGRS2140 

DO  85  1=1, N  DGRS2150 

85  PW(I)  =  Y(I,1)+R* (H*SAVE2 (I) -Y(I,2) )  DGRS2160 

CALL  FCN  (N,T, PW,SAVE1, eprime, eprime2)  + 

NFE  =  NFE+1  DGRS2180 

HL0  =  H*EL(1)  DGRS2190 

DO  90  1=1, N  DGRS2200 

R0  =  H*SAVE2 (I) -Y(I,2)  DGRS2210 

PW(I)  =  1.D0  DGRS2220 

D=  .1D0*R0-H* (SAVEl(I)  -SAVE2 (I)  )                               DGRS2230 

SAVEl(I)  =  0.D0  DGRS2240 

IF  (DABS(RO) .LT.UROUND*YMAX(I) )  GO  TO  90  DGRS2250 

IF  (DABS(D) .EQ.0.D0)  GO  TO  155  DGRS2260 

PW(I)  =  .1D0*R0/D  DGRS2270 

SAVEl(I)  =  PW(I)*R0  DGRS2280 

9  0  CONTINUE  DGRS2290 

GO  TO  135  DGRS2300 

95  IF  (MITER. NE.0)  GOTO  (125,125,105),  MITER  DGRS2310 

C  DGRS2320 

C  IN  THE  CASE  OF  FUNCTIONAL  ITERATION,  DGRS2330 

C  UPDATE  Y  DIRECTLY  FROM  THE  RESULT   DGRS2340 

C  OF  THE  LAST  FCN  CALL.               DGRS2350 

D  =  0.D0  DGRS2360 

DO  100  1=1, N  DGRS2370 

R  =  H*SAVE2  (I)  -Y(I, 2)  DGRS23  80 

D  =  D+( (R-ERROR(I) ) /YMAX(I) ) **2  DGRS2390 

SAVEl(I)  =  Y(I,1)+EL(1)*R  DGRS2400 

100  ERROR (I)  =  R  DGRS2410 

GO  TO  145  DGRS2420 

C  IN  THE  CASE  OF  THE  CHORD  METHOD,      DGRS243  0 

C  COMPUTE  THE  CORRECTOR  ERROR,  F  SUB  DGRS244  0 

C  (M) ,  AND  SOLVE  THE  LINEAR  SYSTEM    DGRS2450 

C  WITH  THAT  AS  RIGHT-HAND  SIDE  AND  P  DGRS2460 
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C  AS  COEFFICIENT  MATRIX,  USING  THE  DGRS24  7  0 

C  LU  DECOMPOSITION  IF  MITER  =  1  OR  DGRS2480 

C  2.  IF  MITER  =  3,  THE  COEFFICIENT  DGRS2490 

C  H*EL(1)  IN  P  IS  UPDATED.  -}GRS2500 

105  PHLO  =  HLO  DGRS2510 

HLO  =  H*EL(1)  DGRS2520 

IF  (HLO. EQ. PHLO)  GO  TO  115  DGRS2530 

R  =  HLO /PHLO  DGRS2  54  0 

DO  110  1=1, N  DGRS2550 

D  =  1.D0-R* (l.D0-l.D0/PW(I) )  DGRS2560 

IF  (DABS(D) .EQ.O.DO)  GOTO  165  DGRS2570 

110  PW(I)  =  1.D0/D  DGRS2580 

115  DO  120  1=1, N  DGRS2590 

120  SAVEl(I)  =  PW(I) * (H*SAVE2 (I) - (Y (I , 2 ) +ERROR (I) ) )  DGRS2600 

GO  TO  13  5  DGRS2  610 

125  DO  130  1=1, N  DGRS2620 

130  SAVEl(I)  =  H*SAVE2 (I) - (Y (I , 2) +ERROR (I) )  DGRS2630 
IF  (NLC  .EQ.  -1)  GO  TO  131  DGRS2  64  0 
NWK  =  (NLC+NUC+1) *N0+1  DGRS2650 
CALL  LEQT1B(PW,N,NLC,NUC,N0,SAVE1,1,N0,2,PW(NWK) , JER)  DGRS2660 
GO  TO  135  DGRS2670 

131  CALL  LUELMF  (PW, SAVE1 , IPIV, N, NO , SAVED  DGRS2680 
135  D  =  0.D0  DGRS2690 

DO  140  1=1, N  DGRS2700 

ERROR (I)  =  ERR0R(I)+SAVE1 (I)  DGRS2710 

D  =  D+(SAVE1(I) /YMAX(I) ) **2  DGRS2720 

140  SAVEl(I)  =  Y(I,1)+EL(1)  *ERROR(I)  DGRS2730 

C  TEST  FOR  CONVERGENCE.  IF  M.GT.0,  THE  DGRS274  0 

C  SQUARE  OF  THE  CONVERGENCE  RATE  DGRS2750 

C  CONSTANT  IS  ESTIMATED  AS  CRATE,  DGRS2760 

C  AND  THIS  IS  USED  IN  THE  TEST.  DGRS2770 

145  IF  (M.NE.0)  CRATE  =  DMAX1 ( . 9D0*CRATE , D/Dl)  DGRS2780 

IF  (  (D*DMIN1  (1  .DO,  2  .D0*CRATE)  )  .LE.BND)  GO  TO  170  DGRS2790 

Dl  =  D  DGRS2800 

M  =  M+l  DGRS2810 

IF  (M.EQ.3)  GO  TO  150  DGRS2820 

CALL  FCN  (N,T, SAVE1 , SAVE2 , eprime , eprime2)  + 

GO  TO  95  DGRS284  0 

C  THE  CORRECTOR  ITERATION  FAILED  TO  DGRS2  850 

C  CONVERGE  IN  3  TRIES.  IF  PARTIALS  DGRS2860 

C  ARE  INVOLVED  BUT  ARE  NOT  UP  TO  DGRS2  870 

C  DATE,  THEY  ARE  REEVALUATED  FOR  THE  DGRS2  880 

C  NEXT  TRY.  OTHERWISE  THE  Y  ARRAY  IS  DGRS2890 

C  RETRACTED  TO  ITS  VALUES  BEFORE  DGRS2900 

C  PREDICTION,  AND  H  IS  REDUCED,  IF  DGRS2910 

C  POSSIBLE.  IF  NOT,  A  NO -CONVERGENCE  DGRS2920 

C  EXIT  IS  TAKEN.  DGRS293  0 

150  NFE  =  NFE+2  DGRS2940 

IF  (IWEVAL.EQ. -1)  GO  TO  165  DGRS2950 

155  T  =  TOLD  DGRS2960 

RMAX  =  2. DO  DGRS2970 

DO  160  J1=1,NQ  DGRS2980 

DO  160  J2=J1,NQ  DGRS2990 

J  =  (NQ+J1) -J2  DGRS3000 

DO  160  1=1, N  DGRS3010 

160  Y(I, J)  =  Y(I,  J)  -Y(I,  J+l)  DGRS3020 

IF  (DABS (H) .LE.HMIN*1.00001D0)  GO  TO  280  DGRS3030 

RH  =  .25D0  DGRS3040 

IREDO  =  1  DGRS3  05  0 

GO  TO  35  DGRS3  060 

165  IWFVAL  =  MITER  DGRS3  070 

GO 'TO  60  DGRS3080 
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c 
c 
c 
c 
c 
c 
c 


c 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


THE  CORRECTOR  HAS  CONVERGED.  IWEVAL 
IS  SET  TO  -1  IF  PARTIAL 
DERIVATIVES  WERE  USED,  TO  SIGNAL 
THAT  THEY  MAY  NEED  UPDATING  ON 
SUBSEQUENT  STEPS.  THE  ERROR  TEST 
IS  MADE  AND  CONTROL  PASSES  TO 
STATEMENT  190  IF  IT  FAILS. 


170 


175 


IF  (MITER. NE.O)  IWEVAL  =  -1 

NFE  =  NFE+M 

D  =  0.D0 

DO  175  1=1, N 

D  =  D+ (ERROR (I) /YMAX(I)  )  **2 

IF  (D.GT.E)  GO  TO  190 


AFTER  A  SUCCESSFUL  STEP,  UPDATE  THE 
Y  ARRAY.  CONSIDER  CHANGING  H  IF 
IDOUB  =  1.  OTHERWISE  DECREASE 
IDOUB  BY  1 .  IF  IDOUB' IS  THEN  1  AND 
NQ  .LT.  MAXDER,  THEN  ERROR  IS 
SAVED  FOR  USE  IN  A  POSSIBLE  ORDER 
INCREASE  ON  THE  NEXT  STEP.  IF  A 
CHANGE  IN  H  IS  CONSIDERED,  AN 
INCREASE  OR  DECREASE  IN  ORDER  BY 
ONE  IS  CONSIDERED  ALSO.  A  CHANGE 
IN  H  IS  MADE  ONLY  IF  IT  IS  BY  A 
FACTOR  OF  AT  LEAST  1.1.  IF  NOT, 
IDOUB  IS  SET  TO  10  TO  PREVENT 
TESTING  FOR  THAT  MANY  STEPS. 


180 


185 


KFLAG  =  0 

IREDO  =  0 

NSTEP  =  NSTEP+1 

HUSED  =  H 

NQUSED  =  NQ 

DO  180  J=1,L 

DO  180  1=1, N 

Y(I, J)  =  Y(I,  J)+EL(J)  *ERROR(I) 

IF  (IDOUB. EQ.l)  GO  TO  200 

IDOUB  =  IDOUB -1 

IF  (IDOUB. GT.l)  GO  TO  290 

IF  (L.EQ.LMAX)  GO  TO  290 

DO  185  1=1, N 

Y(I,LMAX)  =  ERROR  (I) 

GO  TO  29  0 


C 

C 

c 

C 
C 

c 

C 


THE  ERROR  TEST  FAILED.  KFLAG  KEEPS 
TRACK  OF  MULTIPLE  FAILURES . 
RESTORE  T  AND  THE  Y  ARRAY  TO  THEIR 
PREVIOUS  VALUES,  AND  PREPARE  TO 
TRY  THE  STEP  AGAIN.  COMPUTE  THE 
OPTIMUM  STEP  SIZE  FOR  THIS  OR  ONE 
LOWER  ORDER. 


190 


195 


KFLAG  =  KFLAG- 1 
T  =  TOLD 
DO  195  J1=1,NQ 
DO  195  J2=J1,NQ 

J  =  (NQ+J1) -J2 
DO  195  I=1,N 
Y(I,  J)  =  Y(I, J)  -Y(I,J+1) 
RMAX  =  2. DO 

IF  (DABS(H) .LE.HMIN*1.00001D0)  GO  TO  270 
IF  (KFLAG. LE. -3)  GO  TO  260 
IREDO  =  2 
PR3  =  l.D+20 
GO  TO  210 


DGRS3  090 

DGRS3100 

DGRS3110 

DGRS3120 

DGRS3130 

DGRS3140 

DGRS3150 

DGRS3160 

DGRS3170 

DGRS3180 

DGRS3190 

DGRS3200 

DGRS3210 

DGRS3220 

DGRS3230 

DGRS3240 

DGRS3250 

DGRS3260 

DGRS3270 

DGRS3280 

DGRS3290 

DGRS3300 

DGRS3310 

DGRS3320 

DGRS3330 

DGRS3340 

DGRS3350 

DGRS3360 

DGRS3370 

DGRS3380 

DGRS3390 

DGRS3400 

DGRS3410 

DGRS3420 

DGRS3430 

DGRS3440 

DGRS3450 

DGRS3460 

DGRS3470 

DGRS3480 

DGRS3490 

DGRS3500 

DGRS3510 

DGRS3520 

DGRS3530 

DGRS354  0 

DGRS3550 

DGRS3560 

DGRS3570 

DGRS3580 

DGRS359  0 

DGRS3  600 

DGRS3  610 

DGRS3620 

DGRS3  630 

DGRS364  0 

DGRS3  650 

DGRS3  660 

DGRS3670 

DGRS3680 

DGRS3690 

DGRS3  700 
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REGARDLESS  OF  THE  SUCCESS  OR  FAILURE 
OF  THE  STEP,  FACTORS  PR1 ,  PR2 ,  AND 
PR3  ARE  COMPUTED,  BY  WHICH  H  COULD 
BE  DIVIDED  AT  ORDER  NQ  -  1 ,  ORDER 
NQ,  OR  ORDER  NQ  +  1 ,  RESPECTIVELY. 
IN  THE  CASE  OF  FAILURE,  PR3  = 
1.E20  TO  AVOID  AN  ORDER  INCREASE. 
THE  SMALLEST  OF  THESE  IS 
DETERMINED  AND  THE  NEW  ORDER 
CHOSEN  ACCORDINGLY.  IF  THE  ORDER 
IS  TO  BE  INCREASED,  WE  COMPUTE  ONE 
ADDITIONAL  SCALED  DERIVATIVE. 


200  PR3  =  l.D+20 

IF  (L.EQ.LMAX)  GO  TO  210 

Dl  =  0.D0 

DO  205  1=1, N 
205  Dl  =  Dl+( (ERROR(I) -Y(I,LMAX) ) /YMAX(I) ) **2 

ENQ3  =  .5D0/(L+1) 

PR3  =  ( (Dl/EUP) **ENQ3) *1 . 4D0+1 . 4D- 6 
210  ENQ2  =  .5D0/L 

PR2  =  ( (D/E) **ENQ2) *1 . 2D0+1 . 2D- 6 

PR1  =  l.D+2  0 

IF  (NQ.EQ.l)  GO  TO  220 

D  =  0.D0 

DO  215  1=1, N 
215  D  =  D+(Y(I,L) /YMAX(I)  )  **2 

ENQ1  =  .5D0/NQ 

PR1  =  ( (D/EDN) **ENQ1) *1 . 3D0+1 . 3D- 6 
220  IF  (PR2.LE.PR3)  GO  TO  225 

IF  (PR3.LT. PR1)  GO  TO  235 

GO  TO  230 
225  IF  (PR2.GT.PR1)  GO  TO  230 

NEWQ  =  NQ 

RH  =  1.D0/PR2 

GO  TO  250 
23  0  NEWQ  =  NQ-1 

RH  =  1.D0/PR1 

IF  (KFLAG.NE.0.AND.RH.GT.1.D0)  RH  =  1.D0 

GO  TO  250 
235  NEWQ  =  L 

RH  =  1.D0/PR3 

IF  (RH.LT.1.1D0)  GO  TO  245 

DO  240  1=1, N 
240  Y(I,NEWQ+1)  =  ERROR (I) *EL (L) /L 

GO  TO  255 
245  IDOUB  =  10 

GO  TO  29  0 
250  IF  ( (KFLAG.EQ.0) .AND. (RH.LT.1.1D0) )  GO  TO  245 


C 
C 

c 
c 

c 
c 
c 


IF  (NEWQ.EQ.NQ)  GO  TO  35 
255  NQ  =  NEWQ 
L  =  NQ+1 
IRET  =  2 
GO  TO  15 


IF  THERE  IS  A  CHANGE  OF  ORDER,  RESET 
NQ,  L,  AND  THE  COEFFICIENTS.  IN 
ANY  CASE  H  IS  RESET  ACCORDING  TO 
RH  AND  THE  Y  ARRAY  IS  RESCALED . 
THEN  EXIT  FROM  285  IF  THE  STEP  WAS 
OK,  OR  REDO  THE  STEP  OTHERWISE. 


CONTROL  REACHES  THIS  SECTION  IF  3  OR 
MORE  FAILURES  HAVE  OCCURED .  IT  IS 


DGRS3710 

DGRS3720 

DGRS3730 

DGRS3  740 

DGRS3750 

DGRS3760 

DGRS3770 

DGRS3  780 

DGRS3790 

DGRS3800 

DGRS3810 

DGRS3820 

DGRS3830 

DGRS3  84  0 

DGRS3850 

DGRS3860 

DGRS3870 

DGRS3880 

DGRS3  89  0 

DGRS3900 

DGRS3910 

DGRS3920 

DGRS393  0 

DGRS3940 

DGRS3950 

DGRS3960 

DGRS3970 

DGRS39  80 

DGRS3990 

DGRS4  000 

DGRS4010 

DGRS402  0 

DGRS4  03  0 

DGRS4  04  0 

DGRS4  050 

DGRS4  060 

DGRS4070 

DGRS4  080 

DGRS4  090 

DGRS4100 

DGRS4110 

DGRS4120 

DGRS4130 

DGRS414  0 

DGRS4150 

DGRS4160 

DGRS4170 

DGRS4180 

DGRS4190 

DGRS42  00 

DGRS4210 

DGRS4220 

DGRS42  3  0 

DGRS424  0 

DGRS4  250 

DGRS42  60 

DGRS4270 

DGRS4280 

DGRS4290 

DGRS4  3  00 

DGRS4310 

DGRS4320 
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C  ASSUMED  THAT  THE  DERIVATIVES  THAT  DGRS4330 

C  HAVE  ACCUMULATED  IN  THE  Y  ARRAY     DGRS434  0 

C  HAVE  ERRORS  OF  THE  WRONG  ORDER.     DGRS4  350 

C  HENCE  THE  FIRST  DERIVATIVE  IS       DGRS4360 

C  RECOMPUTED,  AND  THE  ORDER  IS  SET    DGRS437  0 

C  TO  1.  THEN  H  IS  REDUCED  BY  A        DGRS4380 

C  FACTOR  OF  10,  AND  THE  STEP  IS      DGRS4390 

C  RETRIED.  AFTER  A  TOTAL  OF  7         DGRS4400 

C  FAILURES,  AN  EXIT  IS  TAKEN  WITH     DGRS4410 

C  KFLAG  =  -2.                         DGRS4420 

260  IF  (KFLAG. EQ. -7)  GO  TO  275  DGRS4430 

RH  =  . 1D0  DGRS4440 

RH  =  DMAX1 (HMIN/DABS (H) , RH)  DGRS4450 

H  =  H*RH  DGRS4460 

CALL  FCN  (N,T, Y, SAVE1 , eprime , eprime2)  + 

NFE  =  NFE+1  DGRS4480 

DO  265  1=1, N  DGRS4490 

265  Y(I, 2)  =  H*SAVE1(I)  DGRS4500 

IWEVAL  =  MITER  DGRS4510 

IDOUB  =  10  DGRS4520 

IF  (NQ.EQ.l)  GO  TO  50  DGRS4530 

NQ  =  1  DGRS454  0 

L  =  2  DGRS4550 

IRET  =  3  DGRS4560 

GO  TO  15  DGRS4570 

C  ALL  RETURNS  ARE  MADE  THROUGH  THIS     DGRS4580 

C  SECTION.  H  IS  SAVED  IN  HOLD  TO     DGRS4590 

C  ALLOW  THE  CALLER  TO  CHANGE  H  ON    DGRS4  600 

C  THE  NEXT  STEP.                       DGRS4  610 

270  KFLAG  =  -1  DGRS4620 

GO  TO  290  DGRS4630 

275  KFLAG  =  -2  DGRS4  640 

GO  TO  290  DGRS4650 

280  KFLAG  =  -3  DGRS4660 

GO  TO  290  DGRS4670 

285  RMAX  =  10. DO  DGRS4680 

290  HOLD  =  H  DGRS4690 

JSTART  =  NQ  DGRS4  700 

C- -Diagnostic  Check  of  first  and  second  derivatives  of  E                + 

if (tcum.eq. told) go  to  310  + 

write  (8,  300)  tcum,  step,  y(  1,1)  ,y(2,l),y(3,l),y(4,l),y(5,l)  + 

300  format (lx,ell.4 , lx, 15, 5 (lx,ell .4) )  + 

write (9,305) step, eprime, eprime2  + 

305  format (lx, 15,2 (lx,e20. 13) )  + 

RETURN  DGRS4  710 

END  DGRS4  720 
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c 

c 

c- 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 

c 


IMSL  ROUTINE  NAME 


DGRCS 


C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 


COMPUTER 
LATEST  REVISION 
PURPOSE 
PRECISION/HARDWARE 

REQD.  IMSL  ROUTINES 
NOTATION 

COPYRIGHT 
WARRANTY 


-  IBM/DOUBLE 

-  JANUARY  1,  1978 

NUCLEUS  CALLED  ONLY  BY  IMSL  SUBROUTINE  DGEAR 

-  SINGLE  AND  DOUBLE /H3 2 

-  SINGLE/H3  6,H4  8,H60 

-  NONE  REQUIRED 

-  INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 

-  1978  BY  IMSL,  INC.  ALL  RIGHTS  RESERVED. 

■  IMSL  WARRANTS  ONLY  THAT  IMSL  TESTING  HAS  BEEN 
APPLIED  TO  THIS  CODE.  NO  OTHER  WARRANTY, 
EXPRESSED  OR  IMPLIED,  IS  APPLICABLE. 


INTEGER 

REAL 

DOUBLE  PRECISION 

INTEGER 

REAL 

DATA 


SUBROUTINE  DGRCS   (METH, NQ, EL, TQ,MAXDER) 

SPECIFICATIONS  FOR  ARGUMENTS 
METH,NQ,MAXDER 
TQ(1) 
EL(1) 

SPECIFICATIONS  FOR  LOCAL  VARIABLES 
K 

PERTST(12,2,3) 
PERTST/1. ,1.,2.,1.,  .3158,  .7407E-1, 

1  .1391E-1,  .2182E-2,  .2945E-3,  .3492E-4, 

2  .3  692E-5,  .3524E-6, 1. ,1.,  .5,  .1667, 

3  .4167E-1,7*1. ,2. ,12. ,24. ,37.89, 

4  53.33,70.08,87.97,106.9,126.7, 

5  14  7.4,168.8,191.0,2.0,4.5,7.333, 

6  10. 42, 13. 7, 7*1., 12. 0,24. 0,37. 89, 

7  53.33,70.08,87.97,106.9,126.7, 

8  147.4,168.8,191.0,1. ,3.0,6.0, 

9  9.167,12.5,8*1./ 

FIRST  EXECUTABLE  STATEMENT 
GO  TO  (5, 10) ,  METH 
5  MAXDER  =  12 

GOTO  (15,20,25,30,35,40,45,50,55,60,65,70),  NQ 
10  MAXDER  =  5 


GOTO  (75,80,85,90,95),  NQ 


THE  FOLLOWING  COEFFICIENTS  SHOULD  BE 
DEFINED  TO  MACHINE  ACCURACY.  FOR  A 
GIVEN  ORDER  NQ,  THEY  CAN  BE 
CALCULATED  BY  USE  OF  THE 
GENERATING  POLYNOMIAL  L (T) ,  WHOSE 
COEFFICIENTS  ARE  EL(I) ..  L (T)  = 
EL(1)  +  EL  (2)  *T  +  ...  + 
EL(NQ+1) *T**NQ.  FOR  THE  IMPLICIT 
ADAMS  METHODS,  L(T)  IS  GIVEN  BY 
DL/DT  =  (T+l)*(T+2)*  ... 
MT+NQ-D/K,  L(-l)  =0,  WHERE  K  = 


DGRC0010 
DGRC0020 
■DGRC0030 
DGRC004  0 
DGRC0050 
DGRC0060 
DGRC0070 
DGRC0080 
DGRC0090 
DGRC0100 
DGRC0110 
DGRC0120 
DGRC013  0 
DGRC0140 
DGRC0150 
DGRC0160 
DGRC017  0 
DGRC0180 
DGRC0190 
DGRC0200 
DGRC0210 
DGRC0220 
DGRC023  0 
DGRC024  0 
DGRC0250 
-DGRC0260 
DGRC0270 
DGRC0280 
DGRC0290 
DGRC03  00 
DGRC0310 
DGRC0320 
DGRC0330 
DGRC034  0 
DGRC0350 
DGRC03  60 
DGRC0370 
DGRC0380 
DGRC0390 
DGRC0400 
DGRC0410 
DGRC0420 
DGRC04  3  0 
DGRC044  0 
DGRC0450 
DGRC0460 
DGRC04  70 
DGRC04  80 
DGRC0490 
DGRC0500 
DGRC0510 
DGRC0520 
DGRC0530 
DGRC054  0 
DGRC0550 
DGRC0560 
DGRC0570 
DGRC0580 
DGRC0590 
DGRC0600 
DGRC0610 
DGRC0620 
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c 
c 
c 
c 
c 
c 
c 
c 
c 


FACTORIAL (NQ-1) .  FOR  THE  GEAR 
METHODS,  L(T)  =  (T+l)*(T+2)*  ... 

* (T+NQ) /K,  WHERE  K  = 
FACTORIAL (NQ) * (1  +  1/2  +  ...  + 

1/NQ) .  THE  ORDER  IN  WHICH  THE 
GROUPS  APPEAR  BELOW  IS. .  IMPLICIT 
ADAMS  METHODS  OF  ORDERS  1  TO  12, 
BACKWARD  DIFFERENTIATION  METHODS 

OF  ORDERS  1  TO  5 . 


15  EL(1) 
GO  TO 

20  EL(1 
EL  (3 
GO  TO 

25  EL(1 
EL  (3 
EL  (4 
GO  TO 

30  EL(1 
EL  (3 
EL  (4 
EL  (5 
GO  TO 

35  EL(1 
EL  (3 
EL  (4 
EL  (5 
EL  (6 
GO  TO 

40  EL(1 
EL  (3 
EL  (4 
EL  (5 
EL  (6 
EL  (7 
GO  TO 

45  EL(1 
EL  (3 
EL  (4 
EL  (5 
EL  (6 
EL  (7 
EL  (8 
GO  TO 

50  EL(1 
EL  (3 
EL  (4 
EL  (5 
EL  (6 
EL  (7 
EL  (8 
EL  (9 
GO  TO 

55  EL(1 
EL  (3 
EL(4 
EL  (5 
EL  (6 
EL  (7 
EL  (8 
EL  (9 
EL(10 


=  1 
100 
=  0. 
=  0. 
100 
=  4  . 
=  0. 
=  1. 
100 
=  0. 
■  9. 
=  3. 
=  4  . 
100 
=  3. 
=  1. 
=  4  . 
=  1  . 
=  8. 
100 
=  3. 
=  1. 
=  0. 
=  1  . 
=  0. 

=  1, 

100 

=  3. 
=  1 
=  7  , 
=  2, 
=  4  . 
=  4 
=  1 
100 
=  3 
=  1 
=  8 
=  3 
=  7 
=  1 
=  7 
=  2 
100 
=  2 
=  1 
=  9 
=  4 
=  1 
=  0 
=  1 
=  1 
2 


0D0 

5D0 
,  5D0 

166666666666667D-01 

75D0 

666666666666667D-01 

375D0 

166666666666667D-01 
333333333333333D-01 
166666666666667D-02 

4  86111111111111D-01 
.041666666666667D0 
. 861111111111111D-01 
. 041666666666667D-01 
.333333333333333D-03 

. 298611111111111D-01 
. 141666666666667D+00 
. 625D+00 

.770833333333333D-01 
.025D+00 
.388888888888889D-03 

. 1559 19 3 12 1693 12D- 01 

.225D+00 

. 51 85 185 185185 19D- 01 

.552083333333333D-01 

. 861111111111111D-02 

.861111111111111D-03 

. 9 84 12 69 84 1269 84D- 04 


04224 
29642 
68518 
35763 
77777 
06481 
93650 
48015 


53703 
85714 
51851 
88888 
77777 
48148 
79365 
87301 


70370D-01 
28571D+00 
85185D-01 
88889D-01 
77778D-02 
14815D-02 
07937D-04 
58730D-05 


)  = 


94 8680004409 171D- 01 
35892857142  8571D+00 
76554232804 232 8D- 01 
171875D-01 
113541666666667D-01 
01875D+00 

934523809523810D-03 
1 1607142 8571429D- 04 
755731922398589D-06 


DGRC0630 

DGRC064  0 

DGRC0650 

DGRC0660 

DGRC0670 

DGRC0680 

DGRC0690 

DGRC0700 

DGRC0710 

DGRC0720 

DGRC0730 

DGRC074  0 

DGRC0750 

DGRC0760 

DGRC0770 

DGRC0780 

DGRC0790 

DGRC0800 

DGRC0810 

DGRC0820 

DGRC0830 

DGRC084  0 

DGRC0850 

DGRC0860 

DGRC0870 

DGRC0880 

DGRC0890 

DGRC0900 

DGRC0910 

DGRC0920 

DGRC0930 

DGRC094  0 

DGRC0950 

DGRC0960 

DGRC0970 

DGRC0980 

DGRC0990 

DGRC1000 

DGRC1010 

DGRC1020 

DGRC103  0 

DGRC104  0 

DGRC1050 

DGRC1060 

DGRC1070 

DGRC1080 

DGRC109  0 

DGRC1100 

DGRC1110 

DGRC112  0 

DGRC1130 

DGRC114  0 

DGRC1150 

DGRC1160 

DGRC1170 

DGRC1180 

DGRC1190 

DGRC1200 

DGRC1210 

DGRC1220 

DGRC123  0 

DGRC124  0 
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60 


65 


70 


75 


80 


85 


90 


95 


GO  TO 

100 

EL(1) 

=  2. 

869  7544  64  2  85714D-01 

EL(3) 

=  1. 

4144  8412  69  84127D+00 

EL  (4) 

=  1. 

0772 156084 65609D+00 

EL(5) 

=  4  . 

985670194003527D-01 

EL  (6) 

=  1. 

484375D-01 

EL  (7) 

=  2. 

9 0605709 8765432D- 02 

EL  (8) 

=  3  . 

720238095238095D-03 

EL  (9) 

=  2. 

996858465608466D-04 

EL(10) 

=  1  . 

3778659 61199295D- 05 

EL(ll) 

=  2 

75573192239 8589D- 07 

GO  TO 

100 

EL(1) 

=  2  . 

801895964439367D-01 

EL(3) 

=  1. 

4  644  8412  69  84127D+00 

EL(4) 

=  1. 

1715145502 64550D+00 

EL  (5) 

=  5 

793581900352734D-01 

EL  (6) 

=  1 

883228615520282D-01 

EL(7) 

=  4 

143036265432099D-02 

EL  (8) 

=  6 

2 1 1144 179 894 1 8 0D- 03 

EL(9) 

=  6 

2520 6679 894 1799D- 04 

EL(10) 

=  4 

04 174 0152 8512 64D- 05 

EL(ll) 

=  1 

.  5156525573 19224D- 06 

EL(12) 

=  2 

.  5052 1083 8544 172D- 08 

GO  TO 

100 

EL(1) 

=  2 

742 6554 003 1599 ID -01 

EL(3) 

=  1 

50993 867243 8672D+00 

EL(4) 

=  1 

2  60271164  021164D+00 

EL(5) 

=  6 

.59234182  09  87654D-01 

EL(6) 

=  2 

.  3 0458002 6455027D- 01 

EL  (7) 

=  5 

5  69724  610523222D-02 

EL  (8) 

=  9 

4  394  8412  6984127D-03 

EL(9) 

=  1 

1192 749 6693 1217D- 03 

EL(10) 

=  9.093915343915344D-05 

EL(ll) 

=  4.822530864197531D-06 

EL(12) 

=  1 

.  503 126503 12 6503D- 07 

EL(13) 

=  2 

.087675698786810D-09 

GO  TO 

100 

EL(1) 

=  1 

.  OD+00 

GO  TO 

100 

EL(1) 

=  6 

.666666666666667D-01 

EL  (3) 

=  3 

.333333333333333D-01 

GO  TO 

100 

EL(1) 

=  5 

. 454 545454 54 5455D- 01 

EL(3) 

=  EL(1) 

EL  (4) 

=  9 

.09  09  09  09  09  09  09  ID -02 

GO  TO 

100 

EL(1) 

a  0 

.48D+00 

EL  (3) 

=  0 

. 7D+00 

EL(4) 

=  0 

.2D+00 

EL  (5) 

=  0 

.  02D+00 

GO  TO 

100 

EL(1) 

=  4 

.379562043795620D-01 

EL  (3) 

=  8 

.211678832116788D-01 

EL(4) 

=  3 

.  102189 78102 189 8D- 01 

EL  (5) 

=  5 

.  4 744 52554 74452 6D- 02 

EL  (6) 

=  3 

.  649  63  503  64  9  6350D-03 

100  DO  105  K=l,3 

TQ(K)  =  PERTST(NQ,METH,K) 
105  CONTINUE 

TQ(4)  =  . 5D0*TQ(2) / (NQ+2) 


DGRC1250 

DGRC1260 

DGRC1270 

DGRC12  80 

DGRC1290 

DGRC1300 

DGRC1310 

DGRC1320 

DGRC1330 

DGRC1340 

DGRC1350 

DGRC13  60 

DGRC13  7  0 

DGRC13  80 

DGRC139  0 

DGRC14  00 

DGRC1410 

DGRC1420 

DGRC1430 

DGRC144  0 

DGRC1450 

DGRC14  60 

DGRC147  0 

DGRC14  80 

DGRC1490 

DGRC1500 

DGRC1510 

DGRC1520 

DGRC153  0 

DGRC1540 

DGRC1550 

DGRC1560 

DGRC1570 

DGRC1580 

DGRC1590 

DGRC1600 

DGRC1610 

DGRC1620 

DGRC1630 

DGRC164  0 

DGRC1650 

DGRC16  60 

DGRC1670 

DGRC1680 

DGRC169  0 

DGRC17  00 

DGRC1710 

DGRC1720 

DGRC173  0 

DGRC1740 

DGRC1750 

DGRC17  60 

DGRC1770 

DGRC1780 

DGRC1790 

DGRC1800 

DGRC1810 

DGRC1820 

DGRC183  0 

DGRC184  0 

DGRC185  0 

DGRC186  0 
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RETURN 
END 


DGRC1870 
DGRC1880 
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c 

c 

c- 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 

c 


IMSL  ROUTINE  NAME 


-  DGRPS 


COMPUTER 
LATEST  REVISION 
PURPOSE 
PRECISION/HARDWARE 

REQD.  IMSL  ROUTINES 
NOTATION 

COPYRIGHT 
WARRANTY 


-  IBM/DOUBLE 

-  NOVEMBER  1,  19  84 

-  NUCLEUS  CALLED  ONLY  BY  IMSL  SUBROUTINE  DGEAR 

-  SINGLE  AND  DOUBLE /H3 2 

-  SINGLE/H3  6,H4  8,H60 

-  LUDATF , LEQT1B , UERTST, UGETIO 

-  INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 

-  19  84  BY  IMSL,  INC.   ALL  RIGHTS  RESERVED. 

-  IMSL  WARRANTS  ONLY  THAT  IMSL  TESTING  HAS  BEEN 

APPLIED  TO  THIS  CODE.   NO  OTHER  WARRANTY, 
EXPRESSED  OR  IMPLIED,  IS  APPLICABLE. 


C 
C 


C 
C 
C 
C 
C 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 


SUBROUTINE  DGRPS  (FCN, FCNJ, Y, NO , CON, MITER, YMAX, SAVE1 , SAVE2 , PW, 

EQUIL,IPIV,IER) 

SPECIFICATIONS  FOR  ARGUMENTS 
NO, MITER, IPIV(l) ,IER 

Y(N0,1)  , CON, YMAX (1)  , SAVE1 (1)  , SAVE  2 (1)  ,PW(1)  , 
EQUIL(l) 

SPECIFICATIONS  FOR  LOCAL  VARIABLES 


INTEGER 

DOUBLE  PRECISION 


INTEGER 


REAL 

DOUBLE  PRECISION 

lr 

COMMON  /DBAND/ 
COMMON  /GEAR/ 


NC ,  MFC ,  KFLAG ,  JSTART ,  NQUSED  ,  NSTEP  ,  NFE  ,  NJE  ,  NPW , 
NSQ , I , Jl , J , NERROR , NSAVE1 , NSAVE2 , NEQUIL , NY , 
IDUMMY (23), NLIM, I I , IJ , LIM1 , LIM2 , NB , NLC , NUC , NWK 

SDUMMY(4) 
T , H , HMIN , HMAX , EPSC , UROUND , EPS J , HUSED , D , RO , YJ , R 
Dl , D2 , WA, DUMMY (4 0 ) 
NLC , NUC 
T , H , HMIN , HMAX , EPSC , UROUND , EPS J , HUSED , DUMMY , 
SDUMMY , NC , MFC , KFLAG , JSTART , NSQ , NQUSED , NSTEP , 
NFE , NJE , NPW , NERROR, NSAVE1 , NSAVE2 , NEQUIL , NY , 
IDUMMY 

THIS  ROUTINE  IS  CALLED  BY  DGRST  TO 
COMPUTE  AND  PROCESS  THE  MATRIX  P  = 
I  -  H*EL(1)*J  ,  WHERE  J  IS  AN 
APPROXIMATION  TO  THE  JACOBIAN.  J 
IS  COMPUTED,  EITHER  BY  THE  USER- 
SUPPLIED  ROUTINE  FCNJ  IF  MITER  = 
1,  OR  BY  FINITE  DIFFERENCING  IF 
MITER  =  2 .  J  IS  STORED  IN  PW  AND 
REPLACED  BY  P,  USING  CON  = 
-H*EL(1) .  THEN  P  IS  SUBJECTED  TO 
LU  DECOMPOSITION  IN  PREPARATION 
FOR  LATER  SOLUTION  OF  LINEAR 
SYSTEMS  WITH  P  AS  COEFFICIENT 
MATRIX.  IN  ADDITION  TO  VARIABLES 
DESCRIBED  PREVIOUSLY, 
COMMUNICATION  WITH  DGRPS  USES  THE 


DGRP0010 
DGRP0020 
-DGRP0030 
DGRP0040 
DGRP0050 
DGRP0060 
DGRP0070 
DGRP0080 
DGRP0090 
DGRP0100 
DGRP0110 
DGRP0120 
DGRP013  0 
DGRP014  0 
DGRP0150 
DGRP0160 
DGRP0170 
DGRP0180 
DGRP0190 
DGRP0200 
DGRP0210 
DGRP0220 
DGRP0230 
DGRP0240 
DGRP0250 
-DGRP0260 
DGRP0270 
DGRP0280 
DGRP0290 
DGRP03  00 
DGRP0310 
DGRP0320 
DGRP033  0 
DGRP0340 
DGRP0350 
DGRP0360 
DGRP0370 
DGRP0380 
DGRP0390 
,DGRP0400 
DGRP0410 
DGRP0420 
DGRP0430 
DGRP0440 
DGRP04  50 
DGRP0460 
DGRP0470 
DGRP0480 
DGRP0490 
DGRP0500 
DGRP0510 
DGRP0520 
DGRP0530 
DGRP0540 
DGRP0550 
DGRP0560 
DGRP0570 
DGRP05  80 
DGRP0590 
DGRP0600 
DGRP0610 
DGRP062  0 
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IF  (NLC.EQ.-l)  GO  TO  45 

NB  =  NLC+NUC+1 

NWK  =  NB*N0+1 

IF  (MITER. EQ. 2)  GO  TO  15 

NLIM  =  NB*N0 
DO  5  1=1, NLIM 

PW(I)  =  0.0D0 
5  CONTINUE 

CALL  FCNJ(NC,T,Y,PW) 
DO  10  1=1, NLIM 

PW(I)  =  PW(I) *CON 
10  CONTINUE 
GO  TO  35 


FOLLOWING  EPSJ  =  DSQRT (UROUND) , 
USED  IN  THE  NUMERICAL  JACOBIAN 
INCREMENTS . 

FIRST  EXECUTABLE  STATEMENT 

BANDED  JACOBIAN  CASE 


MITER  =  1 


MITER  =  2 


15  D  =  0.0D0 

DO  20  1=1, NC 
20  D  =  D+SAVE2 (I) **2 

R0  =  DABS (H) *DSQRT(D) *1.0D+03*UROUND 
DO  3  0  J=1,NC 
YJ  =  Y(J,1) 
R  =  EPSJ*YMAX(J) 
R  =  DMAX1 (R,R0) 
Y(J,1)  =  Y(J,1)+R 
D  =  CON/R 

CALL  FCN(NC,T,Y,  SAVED 
LIM1  =  MAX0 (1, J-NUC) 
LIM2  =  MIN0 (NO, J+NLC) 
DO  25  I=LIM1,LIM2 

I J  =  (J-I+NLC) *N0+I 
PW(IJ)  =  (SAVE1 (I) -SAVE2 (I) ) *D 
25     CONTINUE 

Y(J,1)  =  YJ 

3  0  CONTINUE 

ADD  IDENTITY  MATRIX. 
35  DO  40  1=1, NC 

II  =  NLC*N0+I 

PW(II)  =  PW(II) +1.0D0 

4  0  CONTINUE 

DO  LU  DECOMPOSITION  ON  P 

CALL  LEQT1B(PW,NC,NLC,NUC,N0,EQUIL,1,N0,1,PW(NWK) , IER) 

RETURN 

FULL  JACOBIAN  CASE 
45  IF  (MITER. EQ. 2)  GO  TO  55 

MITER  =  1 

CALL  FCNJ(NC,T,Y,PW) 

DO  50  1=1, NSQ 
50  PW(I)  =  PW(I) *CON 

GO  TO  75 

MITER  =  2 
55  D  =  0.0D0 

DO  60  1=1, NC 
60  D  =  D+SAVE2 (I)**2 

R0  =  DABS (H) *DSQRT(D) *1.0D+03*UROUND 

Jl  =  0 


DGRP0630 

DGRP0640 

DGRP0650 

DGRP0660 

DGRP0670 

DGRP0680 

DGRP0690 

DGRP0700 

DGRP0710 

DGRP0720 

DGRP0730 

DGRP0740 

DGRP0750 

DGRP0760 

DGRP0770 

DGRP0780 

DGRP0790 

DGRP0800 

DGRP0810 

DGRP0820 

DGRP0830 

DGRP0840 

DGRP0850 

DGRP0860 

DGRP0870 

DGRP0880 

DGRP0890 

DGRP0900 

DGRP0910 

DGRP0920 

DGRP0930 

DGRP0940 

DGRP0950 

DGRP0960 

DGRP0970 

DGRP0980 

DGRP0990 

DGRP1000 

DGRP1010 

DGRP1020 

DGRP1030 

DGRP1040 

DGRP105  0 

DGRP1060 

DGRP1070 

DGRP1080 

DGRP1090 

DGRP1100 

DGRP1110 

DGRP1120 

DGRP1130 

DGRP1140 

DGRP1150 

DGRP1160 

DGRP1170 

DGRP1180 

DGRP119  0 

DGRP1200 

DGRP1210 

DGRP1220 

DGRP1230 

DGRP1240 
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DO  70  J=1,NC 
YJ  =  Y(J,1) 
R  =  EPSJ*YMAX(J) 
R  =  DMAX1 (R,R0) 
Y(J,1)  =  Y(J,1)+R 
D  =  CON/R 

CALL  FCN(NC,T,Y, SAVED 
DO  65  1=1, NC 
65     PW(I+J1)  =  (SAVE1 (I) -SAVE2 (I) ) *D 
Y(J,1)  =  YJ 
Jl  =  J1+N0 
7  0  CONTINUE 

ADD  IDENTITY  MATRIX. 
75  J  =  1 

DO  80  1=1, NC 

PW(J)  =  PW(J)+1.0D0 
J  =  J+(N0+1) 
80  CONTINUE 

DO  LU  DECOMPOSITION  ON  P. 

CALL  LUDATF (PW, PW, NC, NO , 0 , Dl , D2 , IPIV, EQUIL, WA, IER) 

RETURN 

END 


DGRP1250 
DGRP12  60 
DGRP1270 
DGRP1280 
DGRP1290 
DGRP1300 
DGRP1310 
DGRP1320 
DGRP1330 
DGRP1340 
DGRP13  50 
DGRP1360 
DGRP1370 
DGRP1380 
DGRP139  0 
DGRP1400 
DGRP1410 
DGRP1420 
DGRP1430 
DGRP1440 
DGRP1450 
DGRP1460 
DGRP1470 
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c 

c 

c- 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 

c 


IMSL  ROUTINE  NAME 


c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


COMPUTER 

LATEST  REVISION 

PURPOSE 

P RE C I S I ON / HARD WARE 

REQD.  IMSL  ROUTINES 
NOTATION 

COPYRIGHT 
WARRANTY 


-  IBM/DOUBLE 

-  JANUARY  1,  1978 
-  NUCLEUS  CALLED  ONLY  BY  IMSL  SUBROUTINE  DGEAR 

-  SINGLE  AND  DOUBLE /H3 2 

-  SINGLE/H3  6,H4  8,H60 

NONE  REQUIRED 

-  INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 

-  1978  BY  IMSL,  INC.  ALL  RIGHTS  RESERVED. 


SUBROUTINE  DGRIN 

INTEGER 

DOUBLE  PRECISION 

INTEGER 


REAL 

DOUBLE  PRECISION 

COMMON  /GEAR/ 


-  DGRIN  DGRI0010 

DGRI0020 

DGRI0030 

DGRI004  0 
DGRI0050 
DGRI0060 
DGRI0070 
DGRI0080 
DGRI0090 
DGRI0100 
DGRI0110 
DGRI0120 
DGRI0130 
DGRI014  0 
DGRI0150 
DGRI0160 
DGRI0170 
DGRI0180 
DGRI0190 
DGRI0200 
DGRI0210 
-  IMSL  WARRANTS  ONLY  THAT  IMSL  TESTING  HAS  BEEN  DGRI0220 
APPLIED  TO  THIS  CODE.  NO  OTHER  WARRANTY,  DGRI0230 
EXPRESSED  OR  IMPLIED,  IS  APPLICABLE.         DGRI024  0 

DGRI0250 

DGRI0260 

DGRI0270 
(TOUT,Y,N0,Y0) 

SPECIFICATIONS  FOR  ARGUMENTS 
NO 
TOUT,Y0 (NO) ,Y(N0,1) 

SPECIFICATIONS  FOR  LOCAL  VARIABLES 
NC , MFC , KFLAG , I , L , J , JSTART , NSQ , NQUSED , NSTEP , 
NFE , NJE , NPW, NERROR, NSAVE1 , NSAVE2 , NEQUIL, NY, 
IDUMMY(23) 
SDUMMY(4) 

T , H , HMIN , HMAX , EPSC , UROUND , EPS J , HUSED , S , SI , 
DUMMY(40) 

T , H , HMIN , HMAX , EPSC , UROUND , EPS J , HUSED , DUMMY , 
SDUMMY , NC , MFC , KFLAG , JSTART , NSQ , NQUSED , NSTEP , 
NFE , NJE , NPW, NERROR, NSAVE1 , NSAVE2 , NEQUIL, NY, 
I DUMMY 

FIRST  EXECUTABLE  STATEMENT 


DO  5  I  =  1,NC 

YO  (I)  =  Y(I,1) 
5  CONTINUE 


DGRI0280 
DGRI0290 
DGRI0300 
DGRI0310 
DGRI0320 
DGRI0330 
DGRI034  0 
DGRI0350 
DGRI0360 
DGRI0370 
DGRI0380 
DGRI0390 
DGRI04  00 
DGRI0410 
DGRI0420 
DGRI0430 
DGRI0440 
DGRI04  50 
DGRI04  60 
THIS  SUBROUTINE  COMPUTES  INTERPOLATEDDGRI04  70 


VALUES  OF  THE  DEPENDENT  VARIABLE 
Y  AND  STORES  THEM  IN  YO .  THE 
INTERPOLATION  IS  TO  THE 
POINT  T  =  TOUT,  AND  USES  THE 
NORDSIECK  HISTORY  ARRAY  Y,  AS 
FOLLOWS . . 

NQ 
Y0(I)   =   SUM   Y(I, J+l) *S**J  , 

J=0 
WHERE  S  =  -(T-TOUD/H. 


L  =  JSTART  +  1 
S  =  (TOUT  -  T) /H 
SI  =  1.0D0 
DO  15  J  =  2,L 
'SI  =  S1*S 


DGRI04  80 
DGRI0490 
DGRI0500 
DGRI0510 
DGRI0520 
DGRI0530 
DGRI054  0 
DGRI0550 
DGRI0560 
DGRI0570 
DGRI0580 
DGRI0590 
DGRI0600 
DGRI0610 
DGRI0620 
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DO  10  I  =  1,NC  DGRI0630 

Y0(I)  =  Y0(I)  +  S1*Y(I,J)                                    DGRI0640 

10     CONTINUE  DGRI0650 

15  CONTINUE  DGRI0660 

RETURN  DGRI0670 

END  DGRI0680 
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IMSL  ROUTINE  NAME 


LUDATF 


COMPUTER 
LATEST  REVISION 
PURPOSE 

USAGE 

ARGUMENTS     A 
LU 


N 
IA 

IDGT 


Dl 


D2 


IPVT 


EQUIL 


WA 


IER 


IBM/DOUBLE 

JANUARY  1,  1978 

L-U  DECOMPOSITION  BY  THE  CROUT  ALGORITHM 
WITH  OPTIONAL  ACCURACY  TEST. 

CALL  LUDATF  (A, LU, N, IA, IDGT, Dl , D2 , IPVT, 
EQUIL, WA, IER) 

INPUT  MATRIX  OF  DIMENSION  N  BY  N  CONTAINING 

THE  MATRIX  TO  BE  DECOMPOSED. 
REAL  OUTPUT  MATRIX  OF  DIMENSION  N  BY  N 
CONTAINING  THE  L-U  DECOMPOSITION  OF  A 
ROWWISE  PERMUTATION  OF  THE  INPUT  MATRIX. 
FOR  A  DESCRIPTION  OF  THE  FORMAT  OF  LU,  SEE 
EXAMPLE . 
INPUT  SCALAR  CONTAINING  THE  ORDER  OF  THE 

MATRIX  A. 
INPUT  SCALAR  CONTAINING  THE  ROW  DIMENSION  OF 
MATRICES  A  AND  LU  EXACTLY  AS  SPECIFIED  IN 
THE  CALLING  PROGRAM. 
INPUT  OPTION. 
IF  IDGT  IS  GREATER  THAN  ZERO,  THE  NON-ZERO 
ELEMENTS  OF  A  ARE  ASSUMED  TO  BE  CORRECT  TO 
IDGT  DECIMAL  PLACES.   LUDATF  PERFORMS  AN 
ACCURACY  TEST  TO  DETERMINE  IF  THE  COMPUTED 
DECOMPOSITION  IS  THE  EXACT  DECOMPOSITION 
OF  A  MATRIX  WHICH  DIFFERS  FROM  THE  GIVEN 
ONE  BY  LESS  THAN  ITS  UNCERTAINTY. 
IF  IDGT  IS  EQUAL  TO  ZERO,  THE  ACCURACY  TEST 

IS  BYPASSED. 
OUTPUT  SCALAR  CONTAINING  ONE  OF  THE  TWO 
COMPONENTS  OF  THE  DETERMINANT.  SEE 
DESCRIPTION  OF  PARAMETER  D2 ,  BELOW. 
OUTPUT  SCALAR  CONTAINING  ONE  OF  THE 

TWO  COMPONENTS  OF  THE  DETERMINANT.  THE 
DETERMINANT  MAY  BE  EVALUATED  AS  (Dl) (2**D2) 
OUTPUT  VECTOR  OF  LENGTH  N  CONTAINING  THE 
PERMUTATION  INDICES .  SEE  DOCUMENT 
(ALGORITHM) . 
OUTPUT  VECTOR  OF  LENGTH  N  CONTAINING 
RECIPROCALS  OF  THE  ABSOLUTE  VALUES  OF 
THE  LARGEST  (IN  ABSOLUTE  VALUE)  ELEMENT 
IN  EACH  ROW. 
ACCURACY  TEST  PARAMETER,  OUTPUT  ONLY  IF 
IDGT  IS  GREATER  THAN  ZERO. 
SEE  ELEMENT  DOCUMENTATION  FOR  DETAILS. 
ERROR  PARAMETER.  (OUTPUT) 
TERMINAL  ERROR 

IER  =  129  INDICATES  THAT  MATRIX  A  IS 
ALGORITHMICALLY  SINGULAR.  (SEE  THE 
CHAPTER  L  PRELUDE) . 
WARNING  ERROR 

IER  =34  INDICATES  THAT  THE  ACCURACY  TEST 
FAILED.   THE  COMPUTED  SOLUTION  MAY  BE  IN 
ERROR  BY  MORE  THAN  CAN  BE  ACCOUNTED  FOR 
BY  THE  UNCERTAINTY  OF  THE  DATA.   THIS 


LUDA0010 
LUDA0020 
-LUDA0030 
LUDA004  0 
LUDA0050 
LUDA0060 
LUDA0070 
LUDA0080 
LUDA0090 
LUDA0100 
LUDA0110 
LUDA0120 
LUDA0130 
LUDA014  0 
LUDA0150 
LUDA0160 
LUDA0170 
LUDA0180 
LUDA0190 
LUDA0200 
LUDA0210 
LUDA0220 
LUDA0230 
LUDA024  0 
LUDA0250 
LUDA0260 
LUDA0270 
LUDA0280 
LUDA0290 
LUDA0300 
LUDA0310 
LUDA0320 
LUDA0330 
LUDA034  0 
LUDA0350 
LUDA0360 
LUDA0370 
LUDA03  80 
LUDA0390 
LUDA04  00 
LUDA0410 
.LUDA04  20 
LUDA043  0 
LUDA044  0 
LUDA04  50 
LUDA04  60 
LUDA04  70 
LUDA04  80 
LUDA0490 
LUDA0500 
LUDA0510 
LUDA052  0 
LUDA053  0 
LUDA054  0 
LUDA0550 
LUDA0560 
LUDA0570 
LUDA0580 
LUDA0590 
LUDA0600 
LUDA0610 
LUDA0620 
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WARNING  CAN  BE  PRODUCED  ONLY  IF  IDGT  IS 
GREATER  THAN  0  ON  INPUT.  SEE  CHAPTER  L 
PRELUDE  FOR  FURTHER  DISCUSSION. 


PRECISION/HARDWARE   -  SINGLE  AND  DOUBLE/H32 

-  SINGLE/H36,H48,H60 

REQD.  IMSL  ROUTINES  -  UERTST, UGETIO 

NOTATION  -  INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 

REMARKS      A  TEST  FOR  SINGULARITY  IS  MADE  AT  TWO  LEVELS: 
1.   A  ROW  OF  THE  ORIGINAL  MATRIX  A  IS  NULL. 
2 .   A  COLUMN  BECOMES  NULL  IN  THE  FACTORIZATION  PROCESS 

COPYRIGHT  -  1978  BY  IMSL,  INC.  ALL  RIGHTS  RESERVED. 

WARRANTY  -  IMSL  WARRANTS  ONLY  THAT  IMSL  TESTING  HAS  BEEN 

APPLIED  TO  THIS  CODE.  NO  OTHER  WARRANTY, 
EXPRESSED  OR  IMPLIED,  IS  APPLICABLE. 


C 
C 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c - -- --- 

c 

SUBROUTINE  LUDATF  (A, LU, N, IA, IDGT, Dl , D2 , IPVT, EQUIL, WA, IER) 
C 

DIMENSION  A(IA,1) ,LU(IA,1) ,IPVT(1) , EQUIL (1) 

DOUBLE  PRECISION  A, LU, Dl , D2 , EQUIL, WA, ZERO, ONE , FOUR, SIXTN, SIXTH, 

*  RN,WREL,BIGA,BIG,P, SUM, AI,WI,T, TEST, Q 
DATA               ZERO , ONE , FOUR, SIXTN, SIXTH/0 . DO , 1 . DO , 4 . DO , 

*  16. DO, .0625D0/ 

C  FIRST  EXECUTABLE  STATEMENT 

C  INITIALIZATION 

IER  =  0 
RN  =  N 
WREL  =  ZERO 
Dl  =  ONE 
D2  =  ZERO 
BIGA  =  ZERO 
DO  10  1=1, N 
BIG  =  ZERO 
DO  5  J=1,N 

P  =  A(I,J) 
LU(I,J)  =  P 
P  =  DABS (P) 

IF  (P  .GT.  BIG)  BIG  =  P 
5     CONTINUE 

IF  (BIG  .GT.  BIGA)  BIGA  =  BIG 
IF  (BIG  .EQ.  ZERO)  GO  TO  110 
EQUIL (I)  =  ONE /BIG 
10  CONTINUE 

DO  105  J=1,N 
JM1  =  J-l 

IF  (JM1  .LT.  1)  GO  TO  40 
C  COMPUTE  U( I, J) ,  1=1,..., J-l 

DO  35  1=1, JM1 

SUM  =  LU(I,J) 
IM1  =1-1 

IF  (IDGT  .EQ.  0)  GO  TO  25 
C  WITH  ACCURACY  TEST 

AI  =  DABS (SUM) 


LUDA0630 

LUDA064  0 

LUDA0650 

LUDA0660 

LUDA0670 

LUDA0680 

LUDA0690 

LUDA0700 

LUDA0710 

LUDA0720 

LUDA07  3  0 

LUDA0740 

LUDA0750 

LUDA0760 

LUDA0770 

, LUDA0780 

LUDA0790 

LUDA0800 

LUDA0810 

LUDA0820 

LUDA0830 

LUDA084  0 

LUDA0850 

-LUDA0860 

LUDA0870 

LUDA0880 

LUDA089  0 

LUDA0900 

LUDA0910 

LUDA0920 

LUDA0930 

LUDA094  0 

LUDA0950 

LUDA0960 

LUDA0970 

LUDA09  80 

LUDA099  0 

LUDA1000 

LUDA1010 

LUDA102  0 

LUDA103  0 

LUDA104  0 

LUDA1050 

LUDA1060 

LUDA1070 

LUDA1080 

LUDA1090 

LUDA1100 

LUDA1110 

LUDA112  0 

LUDA1130 

LUDA114  0 

LUDA1150 

LUDA1160 

LUDA1170 

LUDA1180 

LUDA1190 

LUDA12  00 

LUDA1210 

LUDA1220 

LUDA123  0 

LUDA124  0 


77 


25 


30 

35 

40 


WI  =  ZERO  LUDA1250 

IF  (IM1  .LT.  1)  GO  TO  20  LUDA1260 

DO  15  K=1,IM1  LUDA1270 

T  =  LU(I,K)*LU(K,J)  LUDA1280 

SUM  =  SUM-T  LUDA1290 

WI  =  WI+DABS (T)  LUDA1300 

15        CONTINUE  LUDA1310 

LU(I,J)  =  SUM  LUDA1320 

20        WI  =  WI+DABS  (SUM)  LUDA1330 

IF  (AI  .EQ.  ZERO)  AI  =  BIGA  LUDA134  0 

TEST  =  WI/AI  LUDA1350 

IF  (TEST  .GT.  WREL)  WREL  =  TEST  LUDA1360 

GO  TO  35  LUDA1370 

WITHOUT  ACCURACY  LUDA1380 

IF  (IM1  .LT.  1)  GO  TO  35  LUDA1390 

DO  30  K=1,IM1  LUDA1400 

SUM  =  SUM-LU(I,K)  *LU(K,  J)  LUDA1410 

CONTINUE  LUDA142  0 

LU(I,J)  =  SUM  LUDA1430 

CONTINUE  LUDA1440 

P  =  ZERO  LUDA1450 

COMPUTE  U( J, J)  ANDL(I,J),  I=J+1 , . . . , LUDA1460 

DO  70  I=J,N  LUDA1470 

SUM  =  LU(I,J)  LUDA1480 

IF  (IDGT  .EQ.  0)  GO  TO  55  LUDA1490 

WITH  ACCURACY  TEST  LUDA1500 

AI  =  DABS (SUM)  LUDA1510 

WI  =  ZERO  LUDA1520 

IF  (JM1  .LT.  1)  GO  TO  50  LUDA1530 

DO  45  K=1,JM1  LUDA1540 

T  =  LU(I,K)*LU(K,  J)  LUDA1550 

SUM  =  SUM-T  LUDA1560 

WI  =  WI+DABS (T)  LUDA1570 

45        CONTINUE  LUDA1580 

LU(I,J)  =  SUM  LUDA1590 

50        WI  =  WI+DABS (SUM)  LUDA1600 

IF  (AI  .EQ.  ZERO)  AI  =  BIGA  LUDA1610 

TEST  =  WI/AI  LUDA1620 

IF  (TEST  .GT.  WREL)  WREL  =  TEST  LUDA1630 

GO  TO  65  LUDA164  0 

WITHOUT  ACCURACY  TEST  LUDA1650 

55        IF  (JM1  .LT.  1)  GO  TO  65  LUDA1660 

DO  60  K=1,JM1  LUDA1670 

SUM  =  SUM-LU(I,K) *LU(K, J)  LUDA1680 

60        CONTINUE  LUDA169  0 

LU(I,J)  =  SUM  LUDA1700 

65        Q  =  EQUIL(I)  *DABS  (SUM)  LUDA1710 

IF  (P  .GE.  Q)  GO  TO  70  LUDA1720 

P  =  Q  LUDA1730 

IMAX  =  I  LUDA174  0 

70     CONTINUE  LUDA175  0 

TEST  FOR  ALGORITHMIC  SINGULARITY  LUDA1760 

IF  (RN+P  .EQ.  RN)  GO  TO  110  LUDA1770 

IF  (J  .EQ.  IMAX)  GO  TO  80  LUDA1780 

INTERCHANGE  ROWS  J  AND  IMAX  LUDA1790 

Dl  =  -Dl  LUDA1800 

DO  75  K=1,N  LUDA1810 

P  =  LU(IMAX,K)  LUDA182  0 

LU(IMAX,K)  =  LU(J,K)  LUDA1830 

LU(J,K)  =  P  LUDA1840 

75     CONTINUE  LUDA1850 

'EQUIL(IMAX)  =  EQUIL(J)  LUDA1860 


78 


80 


85 


90 


95 


IPVT(J)  =  IMAX 

Dl  =  D1*LU(J,J) 

IF  (DABS(Dl)  .LE.  ONE)  GO  TO  9  0 

Dl  =  D1*SIXTH 

D2  =  D2+FOUR 

GO  TO  85 

IF  (DABS(Dl)  .GE.  SIXTH)  GO  TO  95 

Dl  =  D1*SIXTN 

D2  =  D2-FOUR 

GO  TO  90 

CONTINUE 

JP1  =  J+l 

IF  (JP1  .GT.  N)  GO  TO  105 


100 

105  CONTINUE 


P  =  LU(J,J) 

DO  100  I=JP1,N 

LU(I,  J)  =  LU(I,  J)  /P 
CONTINUE 


DIVIDE  BY  PIVOT  ELEMENT  U(J,J) 


PERFORM  ACCURACY  TEST 


IF  (IDGT  .EQ.  0)  GO  TO  9005 

P  =  3*N+3 

WA  =  P*WREL 

IF  (WA+10.D0** (-IDGT)  .NE.  WA)  GOTO  9005 

IER  =  34 

GO  TO  9000 


110 


IER  =  129 
Dl  =  ZERO 
D2  =  ZERO 
9  000  CONTINUE 

1 

CALL  UERTST (IER, 6HLUDATF) 
9  005  RETURN 
END 


ALGORITHMIC  SINGULARITY 


PRINT  ERROR 


LUDA1870 
LUDA1880 
LUDA189  0 
LUDA1900 
LUDA1910 
LUDA1920 
LUDA1930 
LUDA194  0 
LUDA1950 
LUDA1960 
LUDA1970 
LUDA19  80 
LUDA199  0 
LUDA2  000 
LUDA2  010 
LUDA202  0 
LUDA2  030 
LUDA2  04  0 
LUDA2050 
LUDA2  060 
LUDA2070 
LUDA2  080 
LUDA2  090 
LUDA2100 
LUDA2110 
LUDA2120 
LUDA213  0 
LUDA214  0 
LUDA2150 
LUDA2160 
LUDA217  0 
LUDA2180 
LUDA2190 
LUDA22  00 
LUDA2210 
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IMSL  ROUTINE  NAME 


LUELMF 


COMPUTER 
LATEST  REVISION 
PURPOSE 

USAGE 
ARGUMENTS     A 


B 
IPVT 


N 
IA 


PRECISION/HARDWARE 

REQD.  IMSL  ROUTINES 
NOTATION 

COPYRIGHT 
WARRANTY 


IBM/DOUBLE 

JANUARY  1,  1978 

ELIMINATION  PART  OF  SOLUTION  OF  AX=B 
(FULL  STORAGE  MODE) 

CALL  LUELMF  (A, B, IPVT, N, IA, X) 

A  =  LU  (THE  RESULT  COMPUTED  IN  THE  IMSL 
ROUTINE  LUDATF)  WHERE  L  IS  A  LOWER 
TRIANGULAR  MATRIX  WITH  ONES  ON  THE  MAIN 
DIAGONAL.  U  IS  UPPER  TRIANGULAR.  L  AND  U 
ARE  STORED  AS  A  SINGLE  MATRIX  A  AND  THE 
UNIT  DIAGONAL  OF  L  IS  NOT  STORED.  (INPUT) 

B  IS  A  VECTOR  OF  LENGTH  N  ON  THE  RIGHT  HAND 
SIDE  OF  THE  EQUATION  AX=B .  (INPUT) 

THE  PERMUTATION  MATRIX  RETURNED  FROM  THE 

IMSL  ROUTINE  LUDATF,  STORED  AS  AN  N  LENGTH 
VECTOR.  (INPUT) 

ORDER  OF  A  AND  NUMBER  OF  ROWS  IN  B.  (INPUT) 

ROW  DIMENSION  OF  A  EXACTLY  AS  SPECIFIED  IN 
THE  DIMENSION  STATEMENT  IN  THE  CALLING 
PROGRAM.  (INPUT) 
THE  RESULT  X.  (OUTPUT) 

-  SINGLE  AND  DOUBLE /H3 2 
SINGLE/H36,H48,H60 

-  NONE  REQUIRED 

•  INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 

-  197  8  BY  IMSL,  INC.  ALL  RIGHTS  RESERVED. 

-  IMSL  WARRANTS  ONLY  THAT  IMSL  TESTING  HAS  BEEN 

APPLIED  TO  THIS  CODE.  NO  OTHER  WARRANTY, 
EXPRESSED  OR  IMPLIED,  IS  APPLICABLE. 


SUBROUTINE  LUELMF  (A, B, IPVT, N, IA,X) 


DIMENSION 
DOUBLE  PRECISION 


DO  5  1=1, N 
5  X(I)  =  B(I) 

IW  =  0 

DO  20  1=1, N 

IP  =  IPVT(I) 
SUM  =  X(IP) 
X(IP)  =  X(I) 

if  (iw  .eq.  o; 

IM1  =  1-1 


A(IA,1)  ,B(1)  ,IPVT(1)  ,X(1) 
A,B,X,SUM 

FIRST  EXECUTABLE  STATEMENT 
SOLVE  LY  =  B  FOR  Y 


GO  TO  15 


LUEF0010 
LUEF0020 
■LUEF0030 
LUEF0040 
LUEF0050 
LUEF0060 
LUEF0070 
LUEF0080 
LUEF0090 
LUEF0100 
LUEF0110 
LUEF0120 
LUEF0130 
LUEF014  0 
LUEF0150 
LUEF0160 
LUEF0170 
LUEF0180 
LUEF0190 
LUEF0200 
LUEF0210 
LUEF0220 
LUEF0230 
LUEF0240 
LUEF0250 
LUEF0260 
LUEF0270 
LUEF0280 
LUEF0290 
LUEF0300 
LUEF0310 
LUEF0320 
LUEF0330 
LUEF0340 
LUEF0350 
LUEF0360 
LUEF0370 
LUEF0380 
LUEF039  0 
LUEF0400 
LUEF0410 
LUEF042  0 
LUEF04  3  0 
LUEF044  0 
LUEF0450 
-LUEF04  60 
LUEF0470 
LUEF04  80 
LUEF0490 
LUEF0500 
LUEF0510 
LUEF0520 
LUEF0530 
LUEF054  0 
LUEF0550 
LUEF0560 
LUEF0570 
LUEF0580 
LUEF0590 
LUEF0600 
LUEF0610 
LUEF0620 


80 


DO  10  J=IW,IM1  LUEF0630 

SUM  =  SUM-ACE,  J)  *X(J)  LUEF0640 

10     CONTINUE  LUEF0650 

GO  TO  20  LUEF0660 

15     IF  (SUM  .NE.  0.D0)  IW  =  I  LUEF0670 

20  X(I)  =  SUM  LUEF0680 

SOLVE  UX  =  Y  FOR  X  LUEF069  0 

DO  30  IB=1,N  LUEF0700 

I  =  N+l-IB  LUEF0710 

IP1  =  1+1  LUEF0720 

SUM  =  X(I)  LUEF073  0 

IF  (IP1  .GT.  N)  GO  TO  30  LUEF0740 

DO  25  J=IP1,N  LUEF0750 

SUM  =  SUM-A(I, J) *X(J)  LUEF0760 

25    CONTINUE  LUEF0770 

30  X(I)  =  SUM/A(I,I)  LUEF0780 

RETURN  LUEF079  0 

END  LUEF0800 
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IMSL  ROUTINE  NAME 


LEQT1B 


COMPUTER 
LATEST  REVISION 
PURPOSE 

USAGE 


ARGUMENTS 


A 

N 

NLC 

NUC 

IA 

B 

M 
IB 

I  JOB 


XL 


IER 


-  IBM/DOUBLE 

-  JANUARY  1,  1978 

-  LINEAR  EQUATION  SOLUTION  -  BAND  STORAGE 

MODE  -  SPACE  ECONOMIZER  SOLUTION 

CALL  LEQT1B  (A, N, NLC, NUC, IA, B, M, IB, IJOB,XL, 
IER) 

-  INPUT/OUTPUT  MATRIX  OF  DIMENSION  N  BY 

(NUC+NLC+1) .  SEE  PARAMETER  I JOB . 
■  ORDER  OF  MATRIX  A  AND  THE  NUMBER  OF  ROWS  IN 
B.  (INPUT) 

-  NUMBER  OF  LOWER  CODIAGONALS  IN  MATRIX  A. 

(INPUT) 

-  NUMBER  OF  UPPER  CODIAGONALS  IN  MATRIX  A. 

(INPUT) 

-  ROW  DIMENSION  OF  MATRIX  A  EXACTLY  AS 

SPECIFIED  IN  THE  DIMENSION  STATEMENT  IN  THE 
CALLING  PROGRAM.  (INPUT) 

-  INPUT/OUTPUT  MATRIX  OF  DIMENSION  N  BY  M. 

ON  INPUT,  B  CONTAINS  THE  M  RIGHT-HAND  SIDES 
OF  THE  EQUATION  AX  =  B .  ON  OUTPUT,  THE 
SOLUTION  MATRIX  X  REPLACES  B.  IF  IJOB  =  1, 
B  IS  NOT  USED. 

-  NUMBER  OF  RIGHT  HAND  SIDES  (COLUMNS  IN  B) . 

(INPUT) 

-  ROW  DIMENSION  OF  MATRIX  B  EXACTLY  AS 

SPECIFIED  IN  THE  DIMENSION  STATEMENT  IN  THE 
CALLING  PROGRAM.  (INPUT) 

-  INPUT  OPTION  PARAMETER.  IJOB  =  I  IMPLIES  WHEN 

1=0,  FACTOR  THE  MATRIX  A  AND  SOLVE  THE 
EQUATION  AX  =  B.  ON  INPUT,  A  CONTAINS  THE 
COEFFICIENT  MATRIX  OF  THE  EQUATION  AX  =  B 
WHERE  A  IS  ASSUMED  TO  BE  AN  N  BY  N  BAND 
MATRIX.  A  IS  STORED  IN  BAND  STORAGE  MODE 
AND  THEREFORE  HAS  DIMENSION  N  BY 
(NLC+NUC+1) .  ON  OUTPUT,  A  IS  REPLACED 
BY  THE  U  MATRIX  OF  THE  L-U  DECOMPOSITION 
OF  A  ROWWISE  PERMUTATION  OF  MATRIX  A.  U 
IS  STORED  IN  BAND  STORAGE  MODE. 

1=1,  FACTOR  THE  MATRIX  A.  A  CONTAINS  THE 
SAME  INPUT/OUTPUT  INFORMATION  AS  IF 
IJOB  =  0 . 

1=2,  SOLVE  THE  EQUATION  AX  =  B .  THIS 
OPTION  IMPLIES  THAT  LEQT1B  HAS  ALREADY 
BEEN  CALLED  USING  IJOB  =  0  OR  1  SO  THAT 
THE  MATRIX  A  HAS  ALREADY  BEEN  FACTORED. 
IN  THIS  CASE,  OUTPUT  MATRICES  A  AND  XL 
MUST  HAVE  BEEN  SAVED  FOR  REUSE  IN  THE 
CALL  TO  LEQT1B. 

-  WORK  AREA  OF  DIMENSION  N*(NLC+1) .  THE  FIRST 

NLC*N  LOCATIONS  OF  XL  CONTAIN  COMPONENTS  OF 
THE  L  MATRIX  OF  THE  L-U  DECOMPOSITION  OF  A 
ROWWISE  PERMUTATION  OF  A.  THE  LAST  N 
LOCATIONS  CONTAIN  THE  PIVOT  INDICES . 

-  ERROR  PARAMETER.  (OUTPUT) 


LE1B0010 
LE1B0020 
-LE1B0030 
LE1B0040 
LE1B0050 
LE1B0060 
LE1B0070 
LE1B0080 
LE1B0090 
LE1B0100 
LE1B0110 
LE1B0120 
LE1B0130 
LE1B0140 
LE1B0150 
LE1B0160 
LE1B0170 
LE1B0180 
LE1B0190 
LE1B0200 
LE1B0210 
LE1B0220 
LE1B0230 
LE1B0240 
LE1B0250 
LE1B0260 
LE1B0270 
LE1B0280 
LE1B0290 
LE1B0300 
LE1B0310 
LE1B0320 
LE1B0330 
LE1B0340 
LE1B0350 
LE1B0360 
LE1B0370 
LE1B0380 
, LE1B0390 
LE1B0400 
LE1B0410 
LE1B0420 
LE1B0430 
LE1B0440 
LE1B0450 
LE1B0460 
LE1B0470 
LE1B0480 
LE1B0490 
LE1B0500 
LE1B0510 
LE1B0520 
LE1B0530 
LE1B0540 
LE1B0550 
LE1B0560 
LE1B0570 
LE1B0580 
LE1B0590 
LE1B0600 
LE1B0610 
LE1B0620 
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TERMINAL  ERROR 

IER  =129  INDICATES  THAT  MATRIX  A  IS 
ALGORITHMICALLY  SINGULAR.  (SEE  THE 
CHAPTER  L  PRELUDE) . 


PRECIS ION /HARDWARE   -  SINGLE  AND  DOUBLE/H32 

-  SINGLE/H3  6,H4  8,H60 

REQD.  IMSL  ROUTINES  -  UERTST, UGETIO 


NOTATION 

COPYRIGHT 
WARRANTY 


INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 

1978  BY  IMSL,  INC.  ALL  RIGHTS  RESERVED. 

IMSL  WARRANTS  ONLY  THAT  IMSL  TESTING  HAS  BEEN 
APPLIED  TO  THIS  CODE.  NO  OTHER  WARRANTY, 
EXPRESSED  OR  IMPLIED,  IS  APPLICABLE. 


SUBROUTINE  LEQT1B  (A, N, NLC , NUC, IA, B, M, IB , IJOB , XL, IER) 


DIMENSION 
DOUBLE  PRECISION 
DATA 

IER  =  0 
JBEG  =  NLC+1 
NLC1  =  JBEG 
IF  (IJOB  .EQ.  2) 
RN  =  N 


1  =  1 

NC  =  JBEG+NUC 
NN  =  NC 
JEND  =  NC 
IF  (N  .EQ 
K  =  1 
P  =  ZERO 
DO  10  J  = 
A(I,K) 


A(IA,1)  ,XL(N,1)  ,B(IB,1) 
A,XL,B,P,Q, ZERO, ONE, RN 
ZERO/0 .DO/, ONE/1 .0D0/ 

FIRST  EXECUTABLE  STATEMENT 


GO  TO  80 


RESTRUCTURE  THE  MATRIX 

FIND  RECIPROCAL  OF  THE  LARGEST 

ABSOLUTE  VALUE  IN  ROW  I 


1  .OR.  NLC  .EQ.  0)  GO  TO  25 


JBEG , JEND 
=  ACE, J) 
Q  =   DABS  (A(I,K)  ) 
IF  (Q  .GT.  P)  P  = 
K  =  K+l 
10  CONTINUE 
IF  (P  .EQ. 
XL(I,NLC1) 
IF  (K  .GT. 
DO  15  J  = 
A(I,J) 
15  CONTINUE 
20  I  =  1+1 

JBEG  =  JBEG-1 

IF  (JEND -JBEG  . EQ .  N)  JEND  =  JEND-1 
IF  (I  .LE.  NLC)  GO  TO  5 
JBEG  =  I 
NN  =  JEND 
25  JEND  =  N-NUC 


ZERO)  GO  TO  13  5 
=  ONE/P 
.  NC)  GO  TO  20 
K,NC 
=  ZERO 


LE1B0630 
LE1B0640 
LE1B0650 
LE1B0660 
LE1B0670 
LE1B0680 
LE1B0690 
LE1B0700 
LE1B0710 
LE1B0720 
LE1B0730 
LE1B0740 
LE1B0750 
LE1B0760 
LE1B0770 
LE1B0780 
LE1B0790 
LE1B0800 
LE1B0810 
LE1B0820 
-LE1B0830 
LE1B0840 
LE1B0850 
LE1B0860 
LE1B0870 
LE1B0880 
LE1B0890 
LE1B0900 
LE1B0910 
LE1B0920 
LE1B0930 
LE1B0940 
LE1B0950 
LE1B0960 
LE1B0970 
LE1B0980 
LE1B0990 
LE1B1000 
LE1B1010 
LE1B1020 
LE1B1030 
LE1B1040 
LE1B1050 
LE1B1060 
LE1B1070 
LE1B1080 
LE1B1090 
LE1B1100 
LE1B1110 
LE1B1120 
LE1B1130 
LE1B1140 
LE1B1150 
LE1B1160 
LE1B1170 
LE1B1180 
LE1B1190 
LE1B1200 
LE1B1210 
LE1B1220 
LE1B1230 
LE1B1240 
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DO  4  0  I  =  JBEG,N 
P  =  ZERO 
DO  30  J  =  1,NN 

Q  =   DABS (A (I, J) ) 
IF  (Q  .GT.  P)  P  =  Q 

3  0     CONTINUE 

IF  (P  .EQ.  ZERO)  GO  TO  135 

XL(I,NLC1)  =  ONE/P 

IF  (I  .EQ.  JEND)  GO  TO  37 

IF  (I  .LT.  JEND)  GO  TO  4  0 

K  =  NN+1 

DO  35  J  =  K,NC 
A (I, J)  =  ZERO 
35     CONTINUE 
37     NN  =  NN-1 

4  0  CONTINUE 

L  =  NLC 

L-U  DECOMPOSITION 
DO  75  K  =  1,N 

P  =   DABS(A(K,1) ) *XL(K,NLC1) 

I  =  K 

IF  (L  .LT.  N)  L  =  L+l 

Kl  =  K+l 

IF  (Kl  .GT.  L)  GO  TO  50 

DO  45  J  =  K1,L 

Q  =  DABS(A(J,1)  )*XL(J,NLC1) 
IF  (Q  .LE.  P)  GO  TO  45 
P  =  Q 
I  =  J 
45     CONTINUE 

50     XL(I,NLC1)  =  XL(K,NLC1) 
XL(K,NLC1)  =  I 

SINGULARITY  FOUND 
Q  =  RN+P 
IF  (Q  .EQ.  RN)  GO  TO  135 

INTERCHANGE  ROWS  I  AND  K 
IF  (K  .EQ.  I)  GO  TO  60 
DO  55  J  =  1,NC 
P  =  A(K,J) 
A(K,J)  =  A(I,J) 
A(I,J)  =  P 
55     CONTINUE 

60     IF  (Kl  .GT.  L)  GO  TO  75 
DO  70  I  =  K1,L 

P  =  A(I,1)  /A(K,1) 
IK  =  I-K 
XL(K1,IK)  =  P 
DO  65  J  =  2,NC 

A(I,J-1)  =  A(I,  J)  -P*A(K,  J) 
65     CONTINUE 

A  (I,  NO  =  ZERO 
70     CONTINUE 
75  CONTINUE 

IF  (I JOB  .EQ.  1)  GO  TO  9  005 

FORWARD  SUBSTITUTION 
80  L  =  NLC 

DO  105  K  =  1,N 

I  =  XL(K,NLC1) 
IF  (I  .EQ.  K)  GO  TO  90 
DO  85  J  =  1,M 
P  =  B(K,  J) 
B(K,J)  =  B(I,J) 


LE1B1250 

LE1B1260 

LE1B1270 

LE1B1280 

LE1B1290 

LE1B1300 

LE1B1310 

LE1B1320 

LE1B1330 

LE1B1340 

LE1B1350 

LE1B1360 

LE1B1370 

LE1B1380 

LE1B1390 

LE1B1400 

LE1B1410 

LE1B1420 

LE1B1430 

LE1B1440 

LE1B1450 

LE1B1460 

LE1B1470 

LE1B1480 

LE1B1490 

LE1B1500 

LE1B1510 

LE1B1520 

LE1B1530 

LE1B1540 

LE1B1550 

LE1B1560 

LE1B1570 

LE1B1580 

LE1B1590 

LE1B1600 

LE1B1610 

LE1B1620 

LE1B1630 

LE1B1640 

LE1B1650 

LE1B1660 

LE1B1670 

LE1B1680 

LE1B1690 

LE1B1700 

LE1B1710 

LE1B1720 

LE1B1730 

LE1B1740 

LE1B1750 

LE1B1760 

LE1B1770 

LE1B1780 

LE1B1790 

LE1B1800 

LE1B1810 

LE1B1820 

LE1B1830 

LE1B184  0 

LE1B1850 

LE1B1860 
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B(I,J)  =  P 
85     CONTINUE 
90     IF  (L  .LT.  N)  L  =  L+l 
Kl  =  K+l 

IF  (Kl  .GT.  L)  GO  TO  105 
DO  100  I  =  K1,L 
IK  =  I-K 
P  =  XL(K1,IK) 
DO  9  5  J  =  1,M 

B(I,  J)  =  B(I,  J)  -P*B(K,  J) 
95        CONTINUE 
100     CONTINUE 
105  CONTINUE 

JBEG  =  NUC+NLC 
DO  125  J  =  1,M 
L  =  1 
Kl  =  N+l 
DO  120  I  =  1,N 
K  =  Kl-I 
P  =  B(K,J) 

IF  (L  .EQ.  1)  GO  TO  115 
DO  110  KK  =  2,L 
IK  =  KK+K 

P  =  P-A(K,KK)  *B(IK-1,  J) 
CONTINUE 
B(K,J)  =  P/A(K,1) 

LE.  JBEG)  L  =  L+l 


BACKWARD  SUBSTITUTION 


110 
115 


IF  (L 
CONTINUE 
125  CONTINUE 

GO  TO  9005 

IER  =  129 

CONTINUE 

CALL  UERTST (IER, 6HLEQT1B) 

RETURN 

END 


120 


135 
9000 

9005 


LE1B1870 

LE1B1880 

LE1B1890 

LE1B1900 

LE1B1910 

LE1B1920 

LE1B1930 

LE1B1940 

LE1B1950 

LE1B1960 

LE1B1970 

LE1B1980 

LE1B1990 

LE1B2000 

LE1B2010 

LE1B2020 

LE1B2030 

LE1B2040 

LE1B2050 

LE1B2060 

LE1B2070 

LE1B2080 

LE1B2090 

LE1B2100 

LE1B2110 

LE1B2120 

LE1B2130 

LE1B2140 

LE1B2150 

LE1B2160 

LE1B2170 

LE1B2180 

LE1B2190 

LE1B2200 

LE1B2210 

LE1B2220 
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IMSL  ROUTINE  NAME 


UERTST 


COMPUTER 

LATEST  REVISION 

PURPOSE 

USAGE 

ARGUMENTS     IER 


NAME 

PRE C I S I ON / HARDWARE 
REQD.  IMSL  ROUTINES 
NOTATION 


■  IBM/SINGLE 

-  JUNE  1,  1982 

PRINT  A  MESSAGE  REFLECTING  AN  ERROR  CONDITION 
CALL  UERTST  (IER, NAME) 

-  ERROR  PARAMETER.  (INPUT) 

IER  =  I+J  WHERE 
I  =  128  IMPLIES  TERMINAL  ERROR  MESSAGE, 
I  =   64  IMPLIES  WARNING  WITH  FIX  MESSAGE, 
I  =   32  IMPLIES  WARNING  MESSAGE. 
J  =  ERROR  CODE  RELEVANT  TO  CALLING 
ROUTINE . 
A  CHARACTER  STRING  OF  LENGTH  SIX  PROVIDING 
THE  NAME  OF  THE  CALLING  ROUTINE.  (INPUT) 

-  SINGLE /ALL 

-  UGETIO,USPKD 

INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 


REMARKS       THE  ERROR  MESSAGE  PRODUCED  BY  UERTST  IS  WRITTEN 
TO  THE  STANDARD  OUTPUT  UNIT.  THE  OUTPUT  UNIT 
NUMBER  CAN  BE  DETERMINED  BY  CALLING  UGETIO  AS 
FOLLOWS . .    CALL  UGETIO ( 1 , NIN, NOUT) . 
THE  OUTPUT  UNIT  NUMBER  CAN  BE  CHANGED  BY  CALLING 
UGETIO  AS  FOLLOWS . . 

NIN  =  0 

NOUT  =  NEW  OUTPUT  UNIT  NUMBER 
CALL  UGETIO (3, NIN, NOUT) 
SEE  THE  UGETIO  DOCUMENT  FOR  MORE  DETAILS. 

COPYRIGHT  -  19  82  BY  IMSL,  INC.  ALL  RIGHTS  RESERVED. 

WARRANTY  -  IMSL  WARRANTS  ONLY  THAT  IMSL  TESTING  HAS  BEEN 

APPLIED  TO  THIS  CODE.  NO  OTHER  WARRANTY, 
EXPRESSED  OR  IMPLIED,  IS  APPLICABLE. 


SUBROUTINE  UERTST  (IER, NAME) 

SPECIFICATIONS  FOR  ARGUMENTS 
IER 
NAME(l) 

SPECIFICATIONS  FOR  LOCAL  VARIABLES 
I , IEQ , IEQDF , IOUNIT , LEVEL , LEVOLD , NAMEQ ( 6 ) , 
NAMSET ( 6 ) , NAMUPK ( 6 ) , NIN, NMTB 
NAMSET/1HU, 1HE, 1HR,1HS, 1HE , 1HT/ 
NAMEQ/6*1H  / 
LEVEL/4/ , IEQDF/0/, IEQ/1H=/ 

UNPACK  NAME  INTO  NAMUPK 
FIRST  EXECUTABLE  STATEMENT 
CALL  USPKD  (NAME, 6, NAMUPK, NMTB) 


INTEGER 
INTEGER 

INTEGER 

r 

DATA 
DATA 
DATA 


UERT0010 
UERT0020 
-UERT0030 
UERT0040 
UERT0050 
UERT0060 
UERT0070 
UERT0080 
UERT0090 
UERT0100 
UERT0110 
UERT0120 
UERT013  0 
UERT0140 
UERT0150 
UERT0160 
UERT0170 
UERT0180 
UERT0190 
UERT0200 
UERT0210 
UERT0220 
UERT0230 
UERT024  0 
UERT0250 
UERT0260 
UERT0270 
UERT0280 
UERT0290 
UERT0300 
UERT0310 
UERT0320 
UERT0330 
UERT034  0 
UERT0350 
UERT0360 
UERT0370 
UERT03  80 
UERT0390 
UERT0400 
UERT0410 
UERT04  20 
UERT0430 
UERT044  0 
UERT0450 
UERT04  60 
UERT0470 
-UERT0480 
UERT0490 
UERT0500 
UERT0510 
UERT0520 
UERT053  0 
UERT054  0 
UERT0550 
UERT0560 
UERT0570 
UERT0580 
UERT0590 
UERT0600 
UERT0610 
UERT0620 
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CALL  UGETI0(1,NIN,I0UNIT) 


GET  OUTPUT  UNIT  NUMBER 


CHECK  IER 


C 

c 
c 
c 


IF  (IER. GT. 999) 

IF  (IER.LT.-32) 

IF  (IER. LE. 128) 

IF  (LEVEL. LT.l) 


GO  TO  25 
GO  TO  55 
GO  TO  5 
GO  TO  30 


CHECK  FOR  UERSET  CALL 


PRINT  TERMINAL  MESSAGE 
IF  (IEQDF.EQ.l)  WRITE (IOUNIT, 35)  IER, NAMEQ, IEQ, NAMUPK 
IF  (IEQDF.EQ.O)  WRITE (IOUNIT, 35)  IER, NAMUPK 
GO  TO  30 
5  IF  (IER.LE.64)  GO  TO  10 
IF  (LEVEL. LT. 2)  GO  TO  30 

PRINT  WARNING  WITH  FIX  MESSAGE 
IF  (IEQDF.EQ.l)  WRITE (IOUNIT, 40)  IER, NAMEQ, IEQ, NAMUPK 
IF  (IEQDF.EQ.O)  WRITE (IOUNIT, 40)  IER, NAMUPK 
GO  TO  30 
10  IF  (IER.LE.32)  GO  TO  15 

PRINT  WARNING  MESSAGE 
IF  (LEVEL. LT. 3)  GO  TO  30 

IF  (IEQDF.EQ.l)  WRITE (IOUNIT, 45)  IER, NAMEQ, IEQ, NAMUPK 
IF  (IEQDF.EQ.O)  WRITE (IOUNIT, 45)  IER, NAMUPK 
GO  TO  3  0 
15  CONTINUE 

DO  20  1=1,6 

IF  (NAMUPK (I) 
20  CONTINUE 

LEVOLD  =  LEVEL 

LEVEL  =  IER 

IER  =  LEVOLD 

IF  (LEVEL. LT.0) 

IF  (LEVEL. GT. 4) 

GO  TO  30 
25  CONTINUE 

IF  (LEVEL. LT. 4) 

IF  (IEQDF.EQ.l) 
IF  (IEQDF.EQ.O) 
30  IEQDF  =  0 

RETURN 

35  FORMAT (19H 

1        20H) 

40  FORMAT (27H 

1        20H) 

45  FORMAT (18H 

1        20H) 

50  FORMAT (20H 


NE.NAMSET(I) )  GO  TO  25 


LEVEL  =  4 
LEVEL  =  4 


GO  TO  30 

PRINT  NON-DEFINED  MESSAGE 
WRITE (IOUNIT, 50)  IER, NAMEQ, IEQ, NAMUPK 
WRITE (IOUNIT, 50)  IER, NAMUPK 


20H) 


***  TERMINAL  ERROR, 1 OX, 7H (IER  =  , 13  , 

FROM  IMSL  ROUTINE  ,6A1,A1,6A1) 
***  WARNING  WITH  FIX  ERROR, 2X, 7H  (IER  = 

FROM  IMSL  ROUTINE  ,6A1,A1,6A1) 
***  WARNING  ERROR, 11X,7H (IER  =  ,  13 , 

FROM  IMSL  ROUTINE  , 6A1 ,  Al , 6A1 ) 
***  UNDEFINED  ERROR, 9X, 7H (IER  =  ,15, 

FROM  IMSL  ROUTINE  ,6A1,A1,6A1) 


13, 


SAVE  P  FOR  P  =  R  CASE 
P  IS  THE  PAGE  NAMUPK 
R  IS  THE  ROUTINE  NAMUPK 


55  IEQDF  =  1 

DO  60  1=1,6 
60  NAMEQ (I)  =  NAMUPK (I) 
65  RETURN 

END 


UERT063  0 

UERT064  0 

UERT0650 

■JERT0660 

UERT067  0 

UERT0680 

UERT069  0 

UERT0700 

UERT0710 

UERT0720 

UERT0730 

UERT074  0 

UERT0750 

UERT0760 

UERT0770 

UERT0780 

UERT079  0 

UERT0800 

UERT0  810 

UERT082  0 

UERT083  0 

UERT084  0 

UERT0850 

UERT0860 

UERT087  0 

UERT0880 

UERT0890 

UERT0900 

UERT0910 

UERT0920 

UERT093  0 

UERT094  0 

UERT0950 

UERT0960 

UERT0970 

UERT09  80 

UERT0990 

UERT1000 

UERT1010 

UERT1020 

UERT1030 

UERT104  0 

UERT1050 

UERT1060 

UERT1070 

UERT1080 

UERT109  0 

UERT1100 

UERT1110 

UERT1120 

UERT1130 

UERT1140 

UERT1150 

UERT1160 

UERT1170 

UERT118  0 

UERT119  0 

UERT1200 
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IMSL  ROUTINE  NAME 


UGETIO 


COMPUTER 
LATEST  REVISION 
PURPOSE 


USAGE 
ARGUMENTS 


IOPT 


NIN 
NOUT 

PRECISION/HARDWARE 


-  IBM/ SINGLE 

-  JUNE  1,  1981 

-  TO  RETRIEVE  CURRENT  VALUES  AND  TO  SET  NEW 

VALUES  FOR  INPUT  AND  OUTPUT  UNIT 
IDENTIFIERS . 

-  CALL  UGETIO (IOPT, NIN, NOUT) 

-  OPTION  PARAMETER.  (INPUT) 

IF  IOPT=l,  THE  CURRENT  INPUT  AND  OUTPUT 
UNIT  IDENTIFIER  VALUES  ARE  RETURNED  IN  NIN 
AND  NOUT,  RESPECTIVELY. 

IF  IOPT=2,  THE  INTERNAL  VALUE  OF  NIN  IS 
RESET  FOR  SUBSEQUENT  USE. 

IF  IOPT=3,  THE  INTERNAL  VALUE  OF  NOUT  IS 
RESET  FOR  SUBSEQUENT  USE. 

-  INPUT  UNIT  IDENTIFIER. 

OUTPUT  IF  IOPT=l,  INPUT  IF  IOPT=2 . 

-  OUTPUT  UNIT  IDENTIFIER. 

OUTPUT  IF  IOPT=l,  INPUT  IF  IOPT=3 . 

-  SINGLE /ALL 


REQD.  IMSL  ROUTINES  -  NONE  REQUIRED 

NOTATION  -  INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 

REMARKS       EACH  IMSL  ROUTINE  THAT  PERFORMS  INPUT  AND /OR  OUTPUT 
OPERATIONS  CALLS  UGETIO  TO  OBTAIN  THE  CURRENT  UNIT 
IDENTIFIER  VALUES.  IF  UGETIO  IS  CALLED  WITH  IOPT=2  OR 
IOPT=3,  NEW  UNIT  IDENTIFIER  VALUES  ARE  ESTABLISHED. 
SUBSEQUENT  INPUT/OUTPUT  IS  PERFORMED  ON  THE  NEW  UNITS. 

COPYRIGHT  -  1978  BY  IMSL,  INC.  ALL  RIGHTS  RESERVED. 

WARRANTY  -  IMSL  WARRANTS  ONLY  THAT  IMSL  TESTING  HAS  BEEN 

APPLIED  TO  THIS  CODE.  NO  OTHER  WARRANTY, 
EXPRESSED  OR  IMPLIED,  IS  APPLICABLE. 


SUBROUTINE  UGETIO (IOPT, NIN, NOUT) 


INTEGER 

INTEGER 
DATA 

IF  (IOPT.EQ.3) 
IF  (IOPT.EQ.2) 
IF  (IOPT.NE.l) 
NIN  =  NIND 
NOUT  =  NOUTD 
GO  TO  9005 


FOR  LOCAL  VARIABLES 


SPECIFICATIONS  FOR  ARGUMENTS 
IOPT, NIN, NOUT 

SPECIFICATIONS 
NIND , NOUTD 
NIND/5/,NOUTD/6/ 

FIRST  EXECUTABLE  STATEMENT 
GO  TO  10 
GO  TO  5 
GO  TO  9005 


UGET0010 
UGET0020 
■UGET0030 
UGET0040 
UGET0050 
UGET0060 
UGET0070 
UGET0080 
UGET0090 
UGET0100 
UGET0110 
UGET0120 
UGET013  0 
UGET0140 
UGET0150 
UGET0160 
UGET0170 
UGET0180 
UGET0190 
UGET0200 
UGET0210 
UGET0220 
UGET0230 
UGET024  0 
UGET0250 
UGET0260 
UGET0270 
UGET0280 
UGET0290 
UGET0300 
UGET0310 
UGET0320 
UGET0330 
UGET034  0 
UGET0350 
UGET0360 
UGET0370 
UGET03  80 
UGET0390 
UGET04  00 
UGET0410 
UGET0420 
UGET0430 
UGET044  0 
UGET04  50 
UGET0460 
UGET0470 
-UGET0480 
UGET0490 
UGET0500 
UGET0510 
UGET052  0 
UGET0530 
UGET054  0 
UGET0550 
UGET0560 
UGET0570 
UGET0580 
UGET0590 
UGET0600 
UGET0610 
UGET062  0 
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5  NIND  =  NIN  UGET063  0 

GO  TO  9005  UGET0640 

10  NOUTD  =  NOUT  UGET0650 

9005  RETURN  UGET0660 

END  UGET0670 
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c 
c 
c- 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c- 

c 

c 


IMSL  ROUTINE  NAME 


-  USPKD 


COMPUTER 
LATEST  REVISION 
PURPOSE 


USAGE 
ARGUMENTS 


-  IBM/SINGLE 

-  NOVEMBER  1,  1984 

-  NUCLEUS  CALLED  BY  IMSL  ROUTINES  THAT  HAVE 

CHARACTER  STRING  ARGUMENTS 

-  CALL  USPKD   ( PACKED, NCHARS , UNPAKD , NCHMTB ) 


PACKED  -  CHARACTER  STRING  TO  BE  UNPACKED. (INPUT) 
NCHARS  -  LENGTH  OF  PACKED.  (INPUT)   SEE  REMARKS. 
UNPAKD  -  INTEGER  ARRAY  TO  RECEIVE  THE  UNPACKED 
REPRESENTATION  OF  THE  STRING.  (OUTPUT) 
NCHMTB  -  NCHARS  MINUS  TRAILING  BLANKS.  (OUTPUT) 


P RE C I S I ON / HARD WARE 
REQD.  IMSL  ROUTINES 
REMARKS   1 


SINGLE /ALL 
NONE 


USPKD  UNPACKS  A  CHARACTER  STRING  INTO  AN  INTEGER  ARRAY 
IN  (Al)  FORMAT. 
2.   UP  TO  129  CHARACTERS  MAY  BE  USED.   ANY  IN  EXCESS  OF 
THAT  ARE  IGNORED. 


COPYRIGHT 
WARRANTY 


-  1984  BY  IMSL,  INC.   ALL  RIGHTS  RESERVED. 

-  IMSL  WARRANTS  ONLY  THAT  IMSL  TESTING  HAS  BEEN 

APPLIED  TO  THIS  CODE.   NO  OTHER  WARRANTY, 
EXPRESSED  OR  IMPLIED,  IS  APPLICABLE. 


SUBROUTINE  USPKD   (PACKED , NCHARS , UNPAKD, NCHMTB) 

SPECIFICATIONS  FOR  ARGUMENTS 
INTEGER  NC , NCHARS , NCHMTB 


LOGICAL*l 
INTEGER*2 


UNPAKD ( 1 ) , PACKED ( 1 ) , LBYTE , LBLANK 
IBYTE , IBLANK 


EQUIVALENCE  ( LBYTE , IBYTE ) 


DATA 
DATA 
DATA 

NCHMTB  =  0 


LBLANK  /1H  / 
IBYTE  /1H  / 
IBLANK  /1H  / 

INITIALIZE  NCHMTB 


IF (NCHARS. LE.O)  RETURN 

NC  =  MINO  (129, NCHARS) 
NWORDS  =  NC*4 
J  =  1 

DO  110  I  =  1, NWORDS, 4 
UNPAKD (I)  =  PACKED (J) 
UNPAKD (1+1)  =  LBLANK 
UNPAKD (1+2)  =  LBLANK 
UNPAKD (1+3)  =  LBLANK 
110  J  =  J+l 


DO' 200  N  =  1, NWORDS, 4 


RETURN  IF  NCHARS  IS  LE  ZERO 


SET  NC=NUMBER  OF  CHARS  TO  BE  DECODED 


CHECK  UNPAKD  ARRAY  AND  SET  NCHMTB 
BASED  ON  TRAILING  BLANKS  FOUND 


USPK0010 
USPK0020 
-USPK0030 
USPK0040 
USPK0050 
USPK0060 
USPK0070 
USPK0080 
USPK0090 
USPK0100 
USPK0110 
USPK0120 
USPK0130 
USPK0140 
USPK0150 
USPK0160 
USPK0170 
USPK0180 
USPK0190 
USPK0200 
USPK0210 
USPK0220 
USPK0230 
USPK0240 
USPK0250 
USPK0260 
USPK0270 
USPK0280 
USPK0290 
USPK0300 
USPK0310 
USPK0320 
USPK0330 
USPK0340 
-USPK0350 
USPK0360 
USPK0370 
USPK0380 
USPK0390 
USPK0400 
USPK0410 
USPK0420 
USPK0430 
USPK0440 
USPK0450 
USPK0460 
USPK0470 
USPK0480 
USPK0490 
USPK0500 
USPK0510 
USPK0520 
USPK0530 
USPK0540 
USPK0550 
USPK0560 
USPK0570 
USPK0580 
USPK0590 
USPK0600 
USPK0610 
USPK0620 
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NN  =  NWORDS  -  N  -  2  USPK0630 

LBYTE  =  UNPAKD(NN)  USPK0640 

IF(IBYTE  .NE.  IBLANK)  GO  TO  210  USPK0650 

200  CONTINUE  USPK0660 

NN  =  0  USPK0670 

210  NCHMTB  =  (NN  +  3)  /  4  USPK0680 

RETURN  USPK0690 

END  USPK0700 
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